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Abstract 

Deriving macroscopic phenomenological laws of irreversible thermodynamics from simple microscopic 
models is one of the tasks of non-equilibrium statistical mechanics. We consider stationary energy 
transport in crystals with reference to simple mathematical models consisting of coupled oscillators on 
a lattice. The role of lattice dimensionality on the breakdown of the Fourier's law is discussed and some 
universal quantitative aspects are emphasized: the divergence of the finite-size thermal conductivity is 
characterized by universal laws in one and two dimensions. Equilibrium and non-equilibrium molecular 



Preprint submitted to Physics Reports 



1 February 2008 



dynamics methods are presented along with a critical survey of previous numerical results. Analytical 
results for the non-equilibrium dynamics can be obtained in the harmonic chain where the role of 
disorder and localization can be also understood. The traditional kinetic approach, based on the 
Boltzmann-Peierls equation is also briefly sketched with reference to one-dimensional chains. Simple 
toy models can be defined in which the conductivity is finite. Anomalous transport in integrable 
nonlinear systems is briefly discussed. Finally, possible future research themes are outlined. 

Key words: Thermal conductivity, classical lattices 
PACS: 63.10.+a, 05.60.-k, 44.10.+i 
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1 Introduction 



The customary macroscopic approach to non-equihbrium phenomena rehes crucially on the 
definition of transport coefficients through phenomenological constitutive equations. Under the 

hypothesis of being close enough to global equilibrium, this is usually accomplished by postulat- 
ing the proportionality among fluxes and thermodynamic forces [1]. For instance when dealing 
with energy transport in a solid one defines the thermal conductivity k through the Fourier's 

Jq = -/«VT , (1) 

where the heat flux Jq is the amount of heat transported through the unit surface per unit time 
and T(x, t) is the local temperature. Such a phenomenological relation was first proposed in 
1808 by J.B.J. Fourier as an attempt to explain the thermal gradient present inside the Earth 
- a problem that had raised a long and controversial debate inside the scientific community 
at that time. Eq. (1) is assumed to be valid close to equilibrium. Actually, the very definition 
of local energy flux Jq(x, and of temperature field ^(x, relies on the local equilibrium 
hypothesis i.e. on the possibility of deflning a local temperature for a macroscopically small but 
microscopically large volume at each location x for each time t. 

The ultimate goal of a complete theory would be to derive an equation like (1) from some 
statistical-mechanics calculation, a task which may be formidably difficult. For insulating crys- 
tals where heat is transported by lattice vibrations, the flrst and most elementary attempt to 
give a microscopic foundation to Fourier's law dates back to Debye. By rephrasing the results of 
the kinetic theory for the (diluted) phonon gas, he found that the thermal conductivity should 
be proportional to Cv£ where C is the heat capacity and v, i are the phonon mean velocity and 
free path, respectively. In 1929, R. Peierls further extended this idea and formulated a Boltz- 
mann equation [2] that shows how anharmonicity is necessary to obtain genuine diffusion of the 
energy through the so-called Umklapp processes. Since then, the Boltzmann- Peierls approach 
became one the cornerstones in the theory of lattice thermal conductivity. Standard methods, 
like the relaxation-time approximation, allow to compute, say, the temperature dependence of 

K. 

From a more fundamental point of view, there are however basic questions that go beyond 
the actual calculation of the transport coefficient. For example, under what condition is local 
equilibrium realized? Can we ensure that a unique non-equilibrium stationary state is attained 
on physically accessible time scales? In this respect, simple mathematical models are an in- 
valuable theoretical playground to provide a more firm foundation to heat conductivity and to 
understand more deeply the hypotheses underlying Eq. (1). Admittedly, this program is still 
nowadays far from being accomplished, at least from a mathematically rigorous point of view 

^ The thermal eonduetivity should be represented in general as a tensor. Here we assume to con- 
sider a simple cubic crystal, in the absence of any external force field. Accordingly, k has a diagonal 
representation with equal diagonal elements. 
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[3]. On the other hand, even in the absence of solvable examples, one can rely on numerical 
simulations as a tool to investigate many of those items. 

As usual in theoretical physics, the guiding criterion of mathematical simplicity leads naturally 
to consider Id or 2d periodic lattices (i.e. chains or planes) of point-like atoms interacting with 
their neighbors through nonlinear forces. The hope is to reproduce realistic thermodynamic 
properties of their three-dimensional counterparts without having to refer to specific struc- 
tures. This brings to the fore the following question: what arc the minimal requirements for a 
dynamical model to satisfy or not Eq. (1)? Although it may appear surprising, this issue has 
been addressed in the literature already in the late 60s without yet receiving a definite answer. 
To a large extent, the present review aims at settling this question by reconciling the very many 
(and sometimes contrasting) numerical results that appeared since then. In fact, several times 
in the past wrong interpretations have been given to the outcome of numerical simulations. In 
the absence of a general theoretical framework, it has been overemphasized the role of deter- 
ministic chaos in ensuring a normal heat transport. Indeed, while ergodicity (implied to some 
extent by a chaotic dynamics) is certainly a necessary condition to establish energy diffusion, 
it has become clear that it is not at all sufficient, as too-rapidly claimed more than a decade 
ago [4]. 

As it is known in the context of fluids [5], much of the difficulties arise from the fact that 
transport coefficients in low spatial dimensions may not exist at all, thus implying a breakdown 
of usual hydrodynamics. In the present context, this usually manifests itself as: (i) a slow decay 
of equilibrium correlations of the heat current; (ii) a divergence of the finite-size conductivity 
k{N, T) in the thermodynamic limit iV — > cxd (where is the number of atoms in the sample). 
One of our concerns will thus be to clarify, through the analysis of several examples, under 
what conditions this should occur. Particular emphasis has been put on the universality of 
quantitative data like the decay law of equilibrium correlation of the flux. 

A by-no-means side issue that is also considered in the present review is the coupling with 
thermal baths. It is only after having properly set the interaction between the system of interest 
and thermal reservoirs that one can be sure that a physically meaningful non-equilibrium regime 
is established. Ideally, a thermal bath is a set of (infinitely) many degrees of freedom, so that it 
can either absorb or release energy, without appreciably changing its own state. Unfortunately, 
such a type of reservoir is very difficult to treat analytically and too much time demanding in 
numerical simulations. Accordingly, various shortcuts have been proposed and are here recalled, 
ranging from stochastic to nonlinear deterministic rules. 

Up to now we have mainly emphasized the theoretical issues that motivate the study of trans- 
port in low- dimensional lattices. Of course, a further relevant motivation is the existence of a 
variety of real systems that could be, at least in principle, effectively described by Id or 2d 
models. For instance, reduced dimensionality has been indeed invoked to explain experiments 
on heat transport in anisotropic crystals [6,7] or magnetic systems [8]. Remarkably, a depen- 
dence of thermal conductivity on the chain length of solid polymers has also been experimen- 
tally observed [9]. More generally, modern experimental techniques [10] allow to directly probe 
the transport properties of semiconductor films [11,12] and single-walled nanotubes [13,14] that 
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markedly display two-dimensional features at low temperatures. Some theoretical investigations 
of thermal conductance for a quantum wire in ballistic [15] and anharmonic [16] regimes have 
been also recently undertaken. Another important example is the problem of heat conductivity 
in quantum spin chains [17]. 



The present work is not simply a review paper in the customary sense: many results were not 
previously published and older results are critically reconsidered. This is of course particularly 
crucial when dealing with numerical data. The plan of the article is the following. In Chapter 
2 we present the simple lattice models that will be considered throughout the review. To be as 
self-contained as possible, we derive the microscopic expression of the heat flux with reference to 
the specific case of one-dimensional systems with nearest-neighbor interactions. The advantage 
of this presentation is twofold: it provides the expression to be referred to in the following and 
allows to understand the hypotheses behind its derivation without having to dwell into more 
involved notations. As already mentioned, an important point for non-equilibrium molecular 
dynamics is the way thermal reservoirs are modeled. Chapter 3 contains a brief survey of some 
simple schemes that have been used in the literature on the topic. The relevant differences 
among the most widely used methods are also discussed. Most of our understanding of energy 
transport in lattices relies on the harmonic approximation for the microscopic dynamics. One 
major advantage of treating the simple harmonic chain is that non-equilibrium properties can 
be derived in a non-perturbative way. This is reviewed in Chapter 4. Harmonic chains are also 
presumably the only class of systems in which the consequences of the presence of quenched 
disorder can be effectively studied. A discussion of chains with isotopic disorder and of the role 
of localization on the heat conduction is given. 



A very sketchy discussion of the two "traditional" approaches, the Boltzmann-Peierls kinetic 
theory for phonons and the Green-Kubo method, is presented in Chapter 5. Since detailed 
accounts exist in many textbooks, we limited ourselves to recall those issues that are relevant 
in the present context. In Chapter 6 we present a detailed account of the many numerical 
studies performed with models where total momentum is conserved. Both equilibrium and 
non-equilibrium simulations are discussed and compared. Chapter 7 is, instead, devoted to 
the "complementary" class of systems where the interaction with an external substrate breaks 
total momentum conservation. This turns out to be a crucial difference in ensuring normal heat 
transport. 



The peculiar behavior of integrable systems is briefly summarized in Chapter 8, while the role 
of dimensionality of the physical space can be appreciated in Chapter 9, where we illustrate 

the behavior of some 2d models. The last chapter is finally devoted to summarizing the key 
points that have been so far understood and, more importantly, to the still open questions. The 
more technical discussions have been confined to the Appendices in order not to downgrade the 
readability of the main body of the text. 
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Fig. 1. A pictorial representation of a chain of A'^ = 4 mutually coupled oscillators in interaction with 
an external substrate and coupled with two thermal reservoirs working at different temperatures. 

2 Definitions 



2.1 Models 



The present Report deals mainly with classical arrays of coupled oscillators. To be more specific, 
we introduce the models for the one- dimensional case. The generalization to two dimension is 
rather straightforward and will be recalled later when needed. 

A schematic setup of the systems that will be mostly studied in the following is drawn in Fig. 1, 
where we have depicted a chain of coupled atoms, the first and the last of which interact 
also with a thermal bath. Let mi and xi be respectively the mass and the position of the l-th 
particle. Only nearest-neighbor interactions will be considered for simplicity. The first class of 
models we wish to consider are defined by an Hamiltonian of the form (pi — mixi) 



N 



pt 

2mi 



+ V{xi+i - xi) 



(2) 



Boundary conditions need also to be specified by defining xo and xn+i- Typical choices are 
periodic, fixed or free boundaries. As only internal forces, that depend on relative positions, 
are present, the total momentum is conserved and thus a zero (Goldstone) mode exist. In 
the harmonic limit, model (2) admits at least a phonon branch whose frequency vanishes for 
vanishing wavenumber. Long-wavelength waves move at the sound velocity and for this reason 
one sometimes refer to (2) as acoustic models. 

An important example is the well-known Lennard- Jones potential 



V{z) = e 



12 



(3) 



In this formulation a is the equilibrium distance and e the well depth. The other example we 
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will often consider is the celebrated Fermi-Pasta-Ulam (FPU) potential [18] 




(4) 



that can be regarded as resulting from an expansion of V close to its equilibrium position z — a. 
Due to its simple algebraic form, the model is computationally very convenient. Two important 
particular cases are worth mentioning: the quadratic plus cubic {g^ = 0) and quadratic plus 
quartic {g^ = 0) potentials that, for historical reasons, are referred to as the FPU-o; and FPU-/3 
models, respectively. In the former one, sufficiently small coupling constant g^ and/or energies 
must be considered to avoid runaway instability of trajectories. 

Models like (2) are a very drastic idealization of a real crystal. Natural low-dimensional lat- 
tice structures are usually embedded in three-dimensional matrices that couple them to the 
environment. Furthermore, artificial arrays of atoms can be constructed by growing them on a 
substrate exerting a pinning force on the atoms in such a way to stabilize the lattice (in gen- 
eral, this is not necessary in the three-dimensional case). At the simplest level of modelization, 
this can be described by adding an external, on-site, potential. For instance, neglecting the 
transverse motion leads to one-dimensional models of the form 



The substrate potential U breaks the invariance xi — > xi + const, of (2) and the total momentum 
is no longer a constant of the motion. Accordingly, all branches of the dispersion relation have 
a gap at zero wavenumber. We therefore refer to (5) as optical models. 

Dimensionless variables will be used throughout whenever possible, especially when reporting 
simulation data. The choice of the most natural units is usually dictated by the particular 
model at hand. For example, for the FPU model with mi = m, it is convenient to set a, m and 
the angular frequency ujq = ^J'gojm to unity. This implies, for instance, that the sound velocity 
OLUo is also unity and that the energy is measured in units of miUQa^. 

2.2 Temperature 

The first problem that has to be solved in order to interpret molecular-dynamics simulations in 
a thermodynamic perspective is the definition of temperature in terms of dynamical variables. 
In appendix A, we show that the problem can be tackled rigorously, although at the expense of 
introducing some technicalities (see also [19,20]). Here below, we follow the traditional approach 
based on the virial theorem. 





T = (u • vn) 



(6) 
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where u is any vector fulfilling the condition V-u = 1, the symbol (■)^ indicates the /i-canonical 
ensemble average, and units are chosen in such a way that the Boltzmann constant is ks — I- 



//-canonical averages are the most appropriate ones whenever an isolated system is numerically 
investigated, but as soon as the system is put in contact with one or more heat baths, canonical 
averages should be considered instead. Fortunately, it is well known, though only partially 
proved, that averages are independent of the chosen ensemble in the thermodynamic limit. 
Nevertheless, when working with finite systems one has to be aware of the existence of finite- 
size corrections as discussed in the celebrated paper by Lebowitz, Percus and Verlet [21]. 

Moreover, in molecular-dynamics simulations, averages are more conveniently computed by 
following single trajectories over time. This is not a problem whenever ergodicity can be invoked, 
so as to ensure that ensemble and time averages are equivalent to one another. However, time 
averages of quantities corresponding to thermodynamic obscrvables have been found to converge 
to the expected ensemble averages even in systems that are known not to be ergodic in a strictly 
mathematical sense (as e.g. the FPU-/3 model at sufficiently small energy values) and even 
when fiuctuations around the mean value are not consistent with equilibrium predictions. This 
suggests that a weaker condition than ergodicity might suffice to ensure the equivalence of time 
and ensemble averages of the physically relevant observables ^ . 

According to Eq. (6), many formally different, but physically equivalent, definitions of tem- 
perature are possible. For instance, the choice u = (0, . . . , 0,pi/N, . . . ,pj^/N) yields the usual 
definition adopted in the canonical ensemble. 



The identification of an optimal definition of temperature to be adopted in numerical studies is 
strictly related to the convergence properties of time-averages. In this sense, it has been observed 
that definitions like the above ones involving only momenta converge always quite rapidly, also 
when the dynamics is weakly chaotic, while definitions involving an explicit dependence on 
space coordinates may converge over much longer time scales [20] . 



^ A possible candidate might be the weak ergodicity criterion, proposed by Khinehin [22], that is 
equivalent to assume a sufficiently fast decay of the time-correlation functions at least for thermody- 
namically relevant observables, like temperature internal energy, specific heat etc. 




(7) 



while the choice u = (0, . . . , 0,pj, 0, . . . , 0) leads to a local definition of temperature. 




(8) 
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2.3 Flux 



The goal of this section is to give a meaningful definition of heat flux [23,24]. This requires some 
care since it involves the transformation of an implicit "mesoscopic" definition into a workable 
microscopic definition. For the sake of simplicity, we discuss the problem with reference to one- 
dimensional systems with nearest-neighbor interactions, the extension to the more general case 
being more technically involved but conceptually equivalent. The heat fiux j{x,t) at time t in 
the spatial position x is nothing but the energy current, implicitly defined by the continuity 
equation 

dh{x,t) ^ dj{x,t) ^ ^ 



where h{x, t) is the energy density. It is important to realize that the energy flux deflned as 

above does not, in general, coincide with heat flux, as the former arises also from macroscopic 
motion [1]. Nonetheless, in solids and one-dimensional fluids no steady motion can occour, so 
that the two fluxes do coincide and we can interchangeably use both names. 

With reference to an ensemble of interacting particles, we can write the microscopic energy 
density as the sum of the isolated contributions located in the instantaneous position of each 
particle 

h{x,t) ^J^^n^i^ - ^n) , (10) 

n 



where d{x) is the Dirac distribution and 



2mr, 



+ U{Xn) + - 



V{Xn+l - Xn) + V{Xn 



X. 



n-l 



is the energy contribution of the nth particle. The flrst two terms on the r.h.s. correspond to the 
kinetic energy and, respectively, to the potential energy f/(a;„) associated with the (possible) 
interaction with an external held. The last term amounts to half of the potential energy of the 
pairwise interactions with the neighboring particles. In a similar way, we can write the heat 
flux as the sum of localized contributions, 

j{x,t) =Y^jJ{x - Xn) . (12) 
n 



The problem amounts therefore to give a deflnition of the local heat flux j„. One should keep 
in mind that the latter quantity has not the same physical dimensions of j{x, t). 

In the limit of small oscillations around the equilibrium position, density fluctuations can be 
neglected and /i„ is equal to the energy density times the lattice spacing a. The time derivative 
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of hr, 



dhn 
dt 



,(13) 



where the prime denotes derivative with respect to the argument and F{x) — —V'{x) can be 
rewritten, with the help of the equations of motion derived from (5) 



(14) 



as 



dhn 
It 



{Xn+l ~l~ Xn)F(^Xn+l 2^n) {Xn ~l~ -^n— l)-^('^n -^n— l) 



(15) 



This equation can, in turn, be rewritten as 

dhji jn jn—1 



dt 



+ 



(16) 



where 



jn = a0n ^a(Xn+l + Xn) F(Xn+l - Xn) 



(17) 



which can thus be interpreted as the local heat flux. 

More in general, if density fluctuations cannot be neglected, a different approach has to be 
followed in order to determine a workable expression for j^- The key idea consists in Fourier 
transforming Eq. (9), 



dh{k,t) 
di 



= -ikj{k,t), 



(18) 



where 



h{k,t) ^ j dxh{x,t)e '''^ j(k,t) ^ j dxj{x,t)e 



■ikx 



(19) 



In fact, according to the idea that the heat flux is a macroscopic observable, one can define it as 
the component of h{k, t) that is proportional to the wave-vector k, i.e. the leading contribution 
over sufficiently large scales. Prom Eq. (10), 



dh{k,t) 
Jt 



/ dhr, 

Vdt 



(20) 
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Prom Eq. (15), the first sum in the r.h.s. can be written by a suitable shift of indexes as 

E = \ - 0n-i)e-'^^" = ^ E ^ne-'"''- (l - e-'=(-''+--")) . (21) 



In the low-A; hmit, the above expression reduces to 



E ^e"''^" ^ -| E(^n+l - Xn)0ne-^'=^" ■ (22) 



n 



By replacing this expression back into Eq. (21) and with the help of Eq. (18), we find that 

jn — '^{^n+l ~ Xn){Xn+l + Xn) Fi^Xn+l — Xn) + Xnhn ■ (23) 

In the limit of small oscillations (compared to the lattice spacing), the second term in the 
above formula can be neglected and Xn — Xn-i — a, so that Eq. (23) reduces to the definition 
(17). Another class of systems for which energy fiuctuations can be neglected is that for which 
the canonical variable "xn" does not describe the longitudinal motion along the chain, but 
corresponds to different degrees of freedom (e.g., transversal oscillations or the rotation of a 
classical spin). In all such cases, the position where the energy /i„ is localized along the chain 
is fixed in time and, accordingly, no term proportional to hn can arise. 

In order to give a flavor of the relative weight of the two contributions to the heat flux, we have 
studied a chain of equal-mass particles interacting throTigh the Lcnnard- Jones potential (3). In 
the low-temperature limit, the chain behaves indeed as a solid, and we expect that the ratio 

ii^n+l ~ ^n)(^n+l + in)F{Xn+l ~ ^n)) 

goes to 0. In the opposite limit T — > oo, the only relevant interaction is the repulsive part of 
the potential that is responsible for elastic collisions, i.e. the system becomes equivalent to a 
hard-point gas. In this limit, the only relevant contribution to the heat flux arises from the 
kinetic term of h„, i.e. 



"ni 

1 



^runxl (25) 



This can be understood by looking at the integral of the flux over the average time between 
consecutive collisions of the same two particles. The flrst term in the r.h.s. of Eq. (23) is not 
negligible only during the inflnitesimal collision time when, however, the distance Xn+i — a;„ 
vanishes. Therefore, in spite of the Dirac S behaviour of the force, its integral contribution 
remains negligible. Analogously, the term XnV{xn+i ~ Xn)-, arising from the potential energy, 
does not contribute, since V remains flnite during the collision. Such a conclusion is confirmed 
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Fig. 2. The ratio rj between the two contributions to the heat flux versus temperature in a chain of 
32 particles with nearest-neighbor Lennard-Jones type interactions. T is the average temperature of 
the two reservoirs while the temperature difference AT is set equal to T/2. 

by numerical simulations which show that the average value of as defined in Eq. (25) coincides 
with the energy flux through the boundaries. 

Altogether, rj is expected to diverge for increasing temperature. This is illustrated in Fig. 2, 
where wc report the data resulting from the simulation of a chain put in contact with two heat 
reservoirs with a temperature difference AT = T/2 and average temperature T. 

Most nonlinear models, aiming at a characterization of solid-like structures have been writ- 
ten with reference to the deviation Qn = Xn — na from the equilibrium position na. This is, 
for instance, the case of the FPU models (4) that can be seen as the result of truncating the 
expansion of the potential energy in powers of (g„ — Qn+i)- The result of introducing the dis- 
placement Qn and getting rid of the actual position Xn, is that physical distances disappear from 
the model and the lattice spacing becomes an arbitrary value. As a consequence, if a is chosen 
too small, one can have unphysical pictures of particles crossing each other and yet keeping 
the same nearest neighbors. A meaningful interpretation of such models is obtained only by 
associating them with a sufficiently large spacing a and thus by computing the heat flux from 
the definition (17). 

Finally, we want to introduce a less symmetric but more compact expression for j„ that can be 
obtained by exploiting the equality (V(g„+i — g„)) = that holds in the stationary regime. By 
determining the derivatives from the equations of motion, we find that 



^ The numerical results have been obtained by implementing a Nose-Hoover thermostat with = 1 
- see the next chapter for further details about the thermostat scheme. 



{Qn+lF{qn+l - Qn)) = (^^-^(^n+l " Qn)) 



(26) 



which allows expressing the average local heat fiow as 




(27) 
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In the stationary regime, it is easily seen that the equahty (^(?n)^) = imphes that 

{qnF{qn+i - qn)) = {qnF{qn - qn-i)) (28) 



The combined use of Eqs. (26,28) for suitable values of n shows that in the bulk 

On) = On-i) := j. (29) 

Accordingly, the average local heat flux is constant along the chain as it should. At the bound- 
aries, one finds that, for whatever choice of the heat bath, the heat flux equals the energy flow 
towards the corresponding reservoir. This can be seen by just writing the balance equation for 
the kinetic energy of the first (last) particle of the chain. 

The quantity that will be mostly used in the following is the total heat flux 

J ^^Jn (30) 
n 



Notice that, from definition (29), (J) = Nj in the stationary regime. From definition (23), J is 
readily recognized to be the one-dimensional version of the general expression (see e.g. [25]) 



(31) 



which is valid for every state of the matter. 

To understand better the physical meaning of the definition given above and for later reference, 
it is useful to consider the case of the simple harmonic chain, i.e. (4) with 5^3 = (74 = and 
ran = iTii- If periodic boundary conditions are assumed, the Hamiltonian is diagonalized by 
passing to the normal coordinates 

1 ^ N N 

Qk-^T.Qi^''^' ' Q-k = Ql, k^-- + l,...,- . (32) 
V-/V ^ ^ 

Considering, for simplicity, definition (17) we get the expression for the total heat flux (30) 

Jh ^ imY^ VkUJkQkQl , (33) 

k 



where — u'^. is the group velocity of phonons and the mode frequencies are given by 



sm 



N 



(34) 
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It is furthermore convenient to introduce the complex amphtudes through the standard 
transformations 

Qk = ^{ake'^^' + a*_,e-'-^'), 4 = '-^ {a,e'^'' - a*_,e-'--' .) . (35) 

A straightforward calculation shows that the heat flux can be expressed as 

Jh = J2 ^kEk , (36) 

k 

where — |ma;^|afep is interpreted as the energy in the A;— th normal mode. This expression, 
originally proposed by Peierls [2], has a simple intuitive interpretation and shows that the 
heat flux is a constant of motion for the harmonic chain at equilibrium. Notice also that its 
equilibrium average {Jh) = as it should: this stems from equipartition of energy (Ek) = ksT 
and from the fact that Vk = —V-k- 



3 Heat baths 

Theoretical and numerical investigations of statistical mechanical systems invariably rely upon 
a suitable modeling of the interaction with thermal reservoirs. At equilibrium, this is usually ac- 
complished by well known methods like micro-canonical molecular dynamics and Monte Carlo 
simulations. Out of equilibrium, the lack of a general theoretical framework forces to define 
meaningful interactions with external thermal baths. The importance of such approach to sim- 
ulations of transport processes in crystalline solids has been already recognized since decades 
[26] 

Conceptually, the correct way to proceed requires considering non-equilibrium states of infinite 
systems. In the context of this review, one could imagine, for example, an infinite chain with 
an initial condition such that all atoms to the left and to the right of some prescribed subset 
(defining the system of interest) are in equilibrium at different temperatures. However, the only 
specific case in which such an approach can be effectively implemented is that of harmonic 
systems [27,28,29,30]. In fact, the degrees of freedom corresponding to the dynamics of the 
reservoirs can be traced out and, as a result, one can prove the existence of stationary non- 
equilibrium states and thereby determine the relevant thermodynamic properties. As soon as 
nonlinear effects are included, it is no longer possible to reduce the evolution of the heat 
baths to a tractable model. One can, nevertheless, study nonlinear chains by assuming that 
nonlinearity is restricted to the system of interest, while still considering linear interactions 
for the semi-infinite chains that correspond to the two reservoirs. By following this approach, 
several results have been obtained (see Ref. [3] and references therein). In particular, it has 
been recently proved the existence of a unique invariant non-equilibrium measure in chains of 
highly nonlinear coupled oscillators subject to arbitrary large temperature gradients [31]. 



14 



If the baths arc modeled with hnear wave equations (the continuum hmit of the harmonic 
lattice), the Hamiltonian dynamics of the full system reduces to the stochastic Markovian 
evolution of a few variables. In the simplest coupling scheme, one can write [3] 

m„qn = -F{qn+l - Qn) + F{qn - Qn-l) + V+Snl + vSnN 

r+ = - A+gi) + ^+ (37) 

r_ = -7-(r_ - X-Qn) + 

where the ^±'s are independent Wiener processes with zero mean and variance 27-i-A-i-T-j-. More- 
over, A± is the coupling strength between the chain and the corresponding bath, while l/7± 
is the relaxation time. This approach is rather recent and we are not aware of any numerical 
study where it is implemented and compared with other methods. 

3.1 Stochastic baths 

A traditional way to implement the interaction with reservoirs amounts to introducing simul- 
taneously random forces and dissipation according to the general prescription of fluctuation- 
dissipation theorem. This could be regarded as the limit case of the previous model when 7-1- 
becomes very large. Consequently, the reservoirs are not affected by the system dynamics. In 
the simple case of an equal-mass chain, this results in the following set of Langevin equations 

mqn = F{qn-q„-i)-F{qn+i-q„) + {^+-X+qn)Sni + {^--X-qn)SnN (38) 

where $,±s are again independent Wiener processes with zero mean and variance 2A±/cs2±- In 
practice, this model too is amenable only to numerical investigation for nonlinear forces. 

Once the nonequilibrium stationary state is attained, one evaluates the average flux j defined 
in(29) and estimates the thermal conductivity as k — \j\L/AT where AT = T+ — T_ . The 
average flux can be obtained also directly from the temperature profile. In fact, direct stochastic 
calculus shows that the average amount of energy exchanged between the first particle and the 
hot reservoir is 

j(A,7V) = ^(T+-rO (39) 
mi 

and an equivalent expression holds on the opposite boundary. This formula states that the heat 
flux is proportional to the difference between the kinetic temperature of the particle in contact 
with the heat bath and the temperature of the reservoir itself. 

A "microscopic" implementation that has been widely used amounts to imagining each reservoir 
as a one-dimensional ideal gas of particles of mass M± interacting with the chain through elastic 
collisions [32] . A simple strategy consists in selecting a random sequence of times ti when each 
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thermostatcd atom collides with a particle of the corresponding reservoir. A natural choice for 
the distribution W{t) of the intervals r between consecutive collisions would be the Poissonian 



iy(r) = iexp(-^), 



(40) 



where r is the average collision time. The kinematics of the collision implies that 



2M. 



{v - qi) 



(41) 



qi + 



m + M. 



+ 



for the left reservoir (an analogous expression holds for the right one). The velocity v of the 
gas particle is a random variable distributed according to the Maxwellian distribution ^ 



In the case M± — m, the procedure reduces to assigning the velocity after the coUision equal 
to random variable v (the particles exchange their velocities). On the other hand, in the limit 
M± <^ m, the interaction with the heat baths becomes Langevin-type as in Eq. (38) with 
A± = 2M±/f . 

This method is computationally very simple, as it avoids the problem of dealing with stochastic 
differential equations: the integration can, in fact, be performed by means of conventional 

algorithms. In particular, the explicit absence of dissipation allows using symplectic integration 
schemes [33,34] between collisions. Furthermore, a physically appealing feature of this approach 
is that damping is not included a priori in the model, but is self-consistently generated by the 
dynamics. 

A related but different scheme consists in determining the collision times from the interaction 
with "thermal walls" suitably placed at the boundaries of the chain. This method has the 
advantage of allowing for the inclusion of pressure effects in the molecular dynamics simulation. 
In this case, the velocity of the thermostated particle is randomized whenever it reaches the 
wall. While the sign of the velocity component normal to the wall has to be chosen in order 
to ensure that it is reflected, its modulus has to be distributed according to a Maxwellian 
distribution at the wall temperature (see Ref . [35] for a discussion of possible pitfalls of different 
choices) . 



^ In one-dimension, it coincides with the Gaussian distribution. 




(42) 
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3.2 Deterministic baths 



In the attempt of providing a self-consistent description of out-of-equilibrium processes, various 
types of deterministic heat baths have been introduced [36]. This was also motivated by the need 
to overcome the difficulties of dealing with stochastic processes. The scheme that has received 
the largest support within molecular-dynamics community is perhaps the so-called Nose-Hoover 
thermostat [37,38]. More precisely, the evolution of the particles in thermal contact with the 
bath a is ruled by the equation 

miin = F{qn - Qn-i) - F{qn+i - gn) - I ^^^'^ if n £ 



where (± are two auxiliary variables modeling the microscopic action of the thermostat, and 
S± denote two sets of N± particles (at the beginning and the end of the chain, respectively) in 
contact with the reservoirs. 

The dynamics of (± is governed by the equation 

.± = 7^ I T.ml-l\ (44) 




where Q± are the thermostat response times. The action of the thermostat can be understood 
in the following terms. Whenever the (kinetic) temperature of the particles in S± is, say, larger 
than T±, (± increases becoming eventually positive. Accordingly, the auxiliary variable acts as 
a dissipation in Eq. (43). Since the opposite occurs whenever the temperature falls below T±, 
this represents altogether a stabilizing feedback around the prescribed temperature. Actually, 
the justification of this scheme rests on a more solid basis. In fact, in Ref. [37,38], it has been 
shown to reproduce the canonical equilibrium distribution. 

In the limit case G — 0, the model reduces to the so-called isokinetic (or Gaussian) thermo- 
stat: the kinetic energy is exactly conserved and the action of the thermal bath is properly 
described without the need to introduce a further dynamical variable, since C± becomes an 
expficit function of the ^'s [36] : 

^ E„GS± Qn [F{qn - qn-l) - F{qn+l - Qn)] , . 

C± -2 ■ (45) 



This latter thermostat scheme can be derived by variational methods after including the non- 
holonomic constraint due to the imposed kinetic energy conservations [36]. 

More generally, it has been shown in Ref. [37,38] that the dynamical equations of this entire 
class of deterministic thermostats possess a Hamiltonian structure in a suitably enlarged phase- 
space. An interesting property that is preserved by the projection onto the usual phase space 
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is time- reversibility. In fact, a simple inspection of the equations reveals that they are invariant 
under time reversal composed with the involution / 

Qn^ -Qn n = 1, . . . , TV , C± ^ -C±. (46) 



This property represents the main reason for the success of this class of thermostats, since 
dissipation is not included a priori, but it rather follows self-consistently from the dynamical 
evolution. In particular, at equilibrium {(±) = 0, indicating that the action of the bath does 
not break microscopic reversibility, while out of equilibrium -\-langleC,-) > and this value 
has been connected with entropy production [39]. 

Another procedure can be defined by combining the idea of a thermostat acting through col- 
lisions at the boundaries with that of a deterministic and reversible rule for the collisions 
themselves. The first context where a scheme of this type has been successfully introduced is 
that of sheared fluids [40]; more recently, van Beijeren [41] has proposed an implementation that 
is suitable for one dimensional systems. The idea is very similar to that of the above discussed 
thermal wall with the crucial difference that the velocity v' after the collision is a deterministic 
and reversible function v' = $(f ) of the initial velocity v. A class of reversible transformations 
is that defined as $ = GRG^^, with R = R^^. A further constraint is that the equilibrium dis- 
tribution be left invariant under the transformation $. The choice adopted in Ref. [90] consists 
in fixing G{v) = exp{—mv^/kBT) and R{x) = 1—x. According to this choice, G transforms the 
Maxwellian distribution into a uniform one on the unit interval, which, in turn, is left invariant 
by R. Although this and the previous choices of reversible thermostats do not correspond to any 
physical mechanism, they offer a convenient framework for the application of dynamical-system 
concepts [39]. 

In the above methods, the system is driven out of equilibrium by a suitable forcing at the 
boundaries: an approach that is aimed at closely reproducing experimental conditions. A dif- 
ferent philosophy [36] consists in introducing an external field acting in the bulk of the system 
to keep it steadily out of equilibrium. The main advantage of this approach is the possibility 
to work with homogeneous systems with, e.g. a uniform temperature along the whole sample. 
In particular, periodic boundary conditions can be enforced thus further reducing finite-size 
effects. This is sometimes referred to as the Evans heat-fiow algorithm and has been applied to 
heat conduction in one-dimensional lattices [42,43,44]. 

More precisely, a fictitious heat field FgDn is added to the equation of motion for the nth 
particle. The coupling Dn must satisfy two conditions: (i) the energy dissipation is proportional 
to jFe, i.e., H — NjFe] (ii) phase-space flux remains divergence-free (this is referred to as the 
adiabatic incompressibility of phase space [36]). Finally, in order to stabilize the dynamics at a 
prescribed temperature, a thermostat rule has to be applied. The resulting equation of motion 
reads as 

miin = F{qn - Qn-l) - F{qn+i - qn) + FeDn - CQu (47) 
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A definition of D„ satisfying the above requirements is [44] 

1 1 ^ 



Moreover, we can consider a Gaussian thermostat: 



1 ^ 



E [^(?n - In-l) - F{<ln+1 ' Qn) + F,Dj] (49) 



Notice that at variance with the previous schemes, the thermostat acts on all particles (this is 
possible since we recall that the temperature is uniform). 

The thermal conductivity can then be determined from the ratio of the heat flux to the applied 
heat field 

K = lim (50) 

where (j) can be equivalently interpreted as a time average or a suitable ensemble average. 
Furthermore, the limit — > is dictated by the need to ensure the validity of the linear 
response theory that is implicitly contained in the initial assumptions. This requirement is all 
the way more substantial in view of the difficulties encountered while working with too large 
heat fields [44] (to some extent this is also true for the Nose-Hoover method descrived above 
[45]). 



3.3 Comparison of different methods 



In all schemes of heat baths there is at least one parameter controUing the coupfing strength: 
let us generically call it g. It can either be the inverse of the average time between subsequent 
collisions, or the dissipation rate A in the Langevin equation, or the inverse of the time-constant 
in the Nose- Hoover scheme. Fixing g is a practical question which is usually solved empirically 
by insuring that different choices do not appreciably affect the outcomes of a simulation. On 
general grounds, one should start by choosing g to be of the order of some typical frequency of 
the system [46]. In the present context, an interesting question immediately arises about the 
dependence of heat transport on g in the various schemes. 

In the case of stochastic reservoirs, the heat fiux vanishes both in the weak- (A — > 0) and 
strong-coupling (A oo) limit. The first implication follows quite easily from the relation (39) 
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and from the observation that (T+ — Ti), cannot increase above (T+ — T_) (more precisely, one 
expects that Ti (T+ + T_)/2 for A ^ 0, since the profile should become increasingly fiat). 
The opposite regime is less trivial: here below, we explain why the same qualitative behavior 
can be found in two different heat-bath schemes. Let us first consider a reservoir acting through 
collisions separated by random intervals r. In this case the coupling constant g is given by the 
inverse of the average collision time r. For small r, one can rely on a pcrturbativc approach 
and write the kinetic energy Ki a time r after the collision as (it is sufficient to consider the 
first particle) 

Ki^{v + Frf (51) 



where v is the random, "initial" velocity. As a result, Ki, on the average, changes by an 
amount that is proportional to r^, since v has zero average. This means that the deviation of 
the (kinetic) temperature from the equilibrium value is also proportional to r^. In order to have 
the energy flux one should multiply this contribution by the number 1/r of collisions per unit 
time. Accordingly, we find that the fiux goes to zero as 1/g in the strong coupling limit. It is 
interesting to notice that the quadratic dependence of the energy could also be inferred from 
the invariance under time reversal of the microscopic equations: after a collision there should 
be no difference between choosing the forward or backward direction for the time axis: as a 
result we have to expect a quadratic behavior! 

Rather similar is the analysis for thermal baths a la Langevin. In that case, for large A, gi can 
be adiabatically eliminated, giving rise to 

£ F 

^i = T + T (52) 



Again, one can see that the average value of the kinetic temperature deviates 0{1/X^) from the 
equilibrium one. 

The simplest way to combine the behavior of the average fiux in the strong and weak coupling 
limit is through the following heuristic formula 



that exhibits the expected A dependence in both limits A ^ and A oo. In Fig. 3 we 
have plotted the fiux for an FPU chain in contact with stochastic Langevin heat baths. Circles 
correspond to direct numerical results, while the solid curve is a fit with formula (53). The 
parameter values turn out to be a = 0.07, bi = 2.45, and 62 = 0.26. The good agreement is to 
be considered as rather incidental, given the heuristic nature of the fitting formula. 

For what concerns the temperature profile, one can see from panel (a) of Fig. 4 that it is 
increasingly fiat for A — > 0. Less obvious is the tendency of the profile to become fiat also in 
the strong coupling- limit. In fact, there is a small but crucial difference for the temperature of 
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Fig. 3. Heat flux dependence on the coupling constant A for an FPU chain in contact with two Langevin 
stochastic reservoirs at temperature T+ = 1.1 and T_ = 0.9. The chain length is N = 128 and fixed 
b.c. have been imposed. 




Fig. 4. Temperature profile in the same set up of the previous figure. In panel (a), the profiles 
correspond to (from bottom to top in the left part) A = 10~^/^, 10~^, 10~^/^, lO'^; panel (b) refers 
to the strong coupling case. Again from bottom to top in the left part, the curves correspond to 
A = 10^4,103/4^10. 

the first and last particle. For A — > the temperature of the extremal particles obviously tends 
to be equal to that in the bulk; on the contrary for A — > oo, such temperatures converge to 
the temperature of the reservoir. In other words the temperature-drop is observed between the 
boundary particle and its neighbor. 

If we repeat the same analysis by using Nose-Hoover thermostats, we find a similar behavior in 
the strong coupling fimit. Indeed, in Fig. 5 we see that the heat fiux vanishes when the response 
time goes to 0. This is the limit of a strictly isokinetic bath and. one is, therefore, bound to 
conclude that such a type of heat baths are unable to sustain heat transport (as long as one 
single particle is thermalized at each extremum). Perhaps more surprising is the opposite limit 
— > oo, since we can see that the heat fiux does not go to zero even though the action of the 
heat baths becomes increasingly slow. 

A partial understanding of this unexpected behavior comes from the observation that the 
variable C, though slowly, reaches the same asymptotic values for arbitrarily large time constants 
0. This is illustrated in Fig. 6 where we have plotted the histogram of ( values for = 10 
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Fig. 5. Heat flux dependence on the coupling constant 1/0 for an FPU chain in contact with two 
Nose-Hoover reservoirs at temperature T+ = 1.1 and r_ = 0.9. The chain length is N = 128 and fixed 
b.c. have been imposed. The line is just a guide for the eyes. 
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Fig. 6. Histogram -P(C) of the left-heat variable in an FPU chain in the same condition as in the 
previous figure and two different values of 9: 10 (solid line), 100 (dashed line) 

(solid line) and 100 (dashed line). Both curves are clean Gaussians centered around the same 
value —0.018.. but with standard deviations differing by a factor 10 (the same factor existing 
between the time constants). In practice, we are led to conclude that, in the limit of — > oo 
the distribution of C values becomes a 5-function centered at a fixed dissipation value that 
depends only on the chain length, the energy and the temperature drop. This result seems to 
be at odds with the fluctuation-dissipation theorem as it seems to suggest that for 0^0 
fluctuations disappear while the dissipation remains flnite. However, one should notice that a 
correct measure of the amount of fluctuations is obtained by integrating over a sufficiently long 
time to let correlations decay: the divergence of the time-constant still ensures the vahdity of 
the fluctuation-dissipation theorem. 

As a result of this analysis we can conclude that values of of order 1 (in the chosen dimen- 
sionless units) are the optimal choices for numerical simulations, since smaller values imply 
smaller heat fluxes, while larger ones would require longer simulation times (in order to ensure 
the decay of correlations) . 
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3-4 Boundary resistances 



Temperature discontinuities usually appear when heat flux is maintained across an interface 
among two substances. This discontinuity is the result of a boundary resistance, generally 
denoted as Kapitza resistance. Its origin is traced back to the "phonon mismatch" between 
two adjacent substances [47]. Such a phenomenon is invariably present in simulations (see 
e.g. Fig. 4 and Fig. 21 below) and may actually reduce the accuracy of the measurements. 
The conductivity evaluated as K^ff = jL/ AT represents an effective transport coefficient that 
includes both boundary and bulk resistances. In practice, one may circumvent the problem 
by referring to a subchain far enough from the ends and compute a bulk conductivity as the 
ratio between J and the actual temperature gradient in the subchain. Notice, however, that 
the advantage of reducing boundary effects may be partly reduced by the increased difficulty 
of dealing with very small temperature differences. 

By denoting with 6T± the temperature jumps at the edges x — and x — L, the imposed 
temperature difference, AT = T+ — T_, can be written as 

L 

AT = 5T+ + 5T_ + J dxVT 



Let us now express the jumps 5T± as [48] 

ST± = ei±VT\^ (54) 



where VT\± is the value of the temperature gradient extrapolated at the boundary {x — 
0, L), £± is the mean free path at the corresponding temperature; the phenomcnological and 
dimensionless constant e measures the coupling strength between the thermostats and the 
system. 

This equation has been numerically tested by plotting 5T/£ versus VT" for different thermostats 
[49,50]. From the slope of the various curves it has been found that in the FPU-/3 model the 
coupling parameter may vary from e = 0.8 for the thermostats adopted in Ref. [50], to e = 2. 
for Nose-Hoover thermostats, to e = 400 for stochastic reservoirs as evidenced in the large 
boundary jumps seen, e.g., in Fig. 4. This reflects the empirical fact that the shape of the 
temperature profile for given T± may depend on the choice of thermostats. Nevertheless, the 
scaling is independent of this choice. 

In the near-equilibrium regime, we can assume ~ 5T_, ~ £_ = £ and the gradient to be 
constant, VT|-|- = VT obtaining 

(2ee + L)VT^ AT. (55) 
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If we further denote by k — j/VT the bulk conductivity we obtain 



3 



(56) 



AT/L 




When the mean free path is much larger than the system size i ^ L (the so-called Casimir 
limit [48]) energy flows almost freely through the system, boundary scattering takes over, the 
system is almost harmonic with Kgff oc L and a flat temperature profile occurs. In the opposite 
case i <^ L, Kgff — n and one is indeed probing the bulk interaction. 



4 Hcirmonic systems 

The simplest and almost unique class of systems for which one can perform analytic calculations 
is represented by harmonic chains. Even though they are characterized by a peculiar dynamics, 
basically because of the integrability of the motion, their behavior can help shedding some 
light on various aspects of heat conductivity. One of the properties of harmonic chains is the 
possibility to decompose the heat flux into the sum of independent contributions associated to 
the various cigcnmodes. This analysis is particularly useful to obtain a deeper insight about 
the role of boundary conditions. 

We first discuss the dynamics of homogeneous chains, since it is possible to obtain an analytic 
expression for the invariant measure in the general case of arbitrary coupling strength. The 
effect of disorder is discussed in the subsequent section, where perturbativc calculations of the 
relevant quantities are illustrated. For completeness, we also recall the localization properties of 
the eigenfunctions and self-averaging properties of several observables such as the temperature 
profile and the heat flux. The approach that we have followed is mainly based on the Fokker- 
Planck equation and simple stochastic calculus. An important alternative approach based on 
transmission coefficients can be found in Ref. [28]. 



4-1 Homogeneous chains 

We consider a homogeneous harmonic chain with fixed boundary conditions in contact with 
stochastic Langevin heat baths. Eq. (38) reduces to 



where we have set A+ = A_ = A to lighten the notations and assumed unitary masses. This 
set of stochastic equations can be solved [51] by passing to a phase-space description, i.e. by 



) 5ni{i+ - Agi) 8nN{i- - Agjv), 



(57) 
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writing down the Liouville equation that, in this case, corresponds to the following Fokker- 
Planck equation 



<9P _ d d^P 



dt ^^^^dxi ' ' 2 dxidx 



(58) 



where here and in the following we adopt the summation convention (i.e., sums over repeated 
indices are understood without explicitly writing down the summation sign); Xi = qi for 1 < 
i < N, Xi = (ji for N < i < 2N. Aij and Dij are elements of the 2N x 2N matrices A and D 
that we write in terms of N x N blocks 



-I 

u'^G AR 



D 



° ^ 

2XkBT{K + r]S) J 



(59) 



where we have introduced the average temperature T = (T+ + T_)/2 and the rescaled temper- 
ature difference rj — (7+ — TJ)/T. Moreover, and I are the null and identity matrices, G is a 
tridiagonal matrix defined as 



Gij — 2Sij — — ^ij+i , 



while R and S are defined as 



Rij — Sij{Sii + Sin), 
Sij — Sij{Sii — Sin) 

The general solution of this equation can be sought of the form 



(60) 



(61) 



where C is the symmetric covariance matrix 

Cij = (xiXj) = J dxP{x) 



(62) 



By replacing the definition of C into Eq. (58), one finds that 

C = D - AC - CA^ 



(63) 



where A''' denotes the transpose of A. Accordingly, the asymptotic stationary solutions can be 
determined from the following equation [52] 

D = AC + CA^ (64) 
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^={^5 I) (65) 



In order to solve the problem, let us write C in terms oi N x N blocks, 

U Z 

.Z^ V 

where the matrices U, V, and Z express the correlations among positions and velocities, 

Uij = {qiqj) , Vij = {qiQj) , Zij = {q^q^) , (66) 

If the temperatures of the two heat baths coincide (i.e. 77 = 0), it can be easily seen that 

Ue = ^G-^ , Ve = A;sri , Ze = (67) 

represent a meaningful solution, since it coincides with the equilibrium Boltzmann distribution 
P{x) ^ ex.p{-H/kBT). 

The derivation of the stationary solution in the out-of-equilibrium case is reported in App. B. 
All relevant correlations can be expressed in terms of the function (p{j) (see Eq. (B.13)) that 
decays exponentially with the rate a defined in Eq. (B.12). As it can be seen by direct inspection 
of correlations, a measures the length over which the boundary reservoirs significantly affect 
the chain dynamics. As expected, a diverges in the weak coupling limit {u — uP' j }? 00) . 

From Eqs. (B.8,B.10), it follows that position-position and velocity- velocity correlations are 
equal for all pairs of particles (i,j) such that i -|- j is constant. The qualitative explanation of 
this property relies on the exponential decay of the boundary effects. In fact, the amplitude 
of, e.g., (gjgj) decreasing exponentially with both i and j has to depend on i + j. Less obvious 
is the left-right antisymmetry {V^ = VN^i^i^N^j^i), which implies that the boundary effects 
are exactly the same for the two thermostats, whatever is their temperature. This is nicely 
reproduced by the temperature profile 



T{z)=T{l + r]Vu) = { 



T+ - yriT(p{l) i = 1 

T[l-W(2i-1)] l<i<N/2 

T[l + r]V(t){2{N - i) - l)\ N/2<i<N ^ ' 

{ T_ + uriT(f){l) i = N 



which exhibits a further unexpected property (see Fig. 7): the temperature is higher in the 
vicinity of the coldest reservoir (the only exception being represented by the first and last par- 
ticles)!. Because of the exponential decay of 0(i), in the bulk, the temperature profile is constant 
as if the system were at equilibrium at temperature T. However, this is only superficially true, 
as position-velocity correlations significantly differ from the equilibrium ones. 

Also the average stationary local flux can be computed explicitly: 

ji^u Zi^i^i = 0(1) (69) 
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Fig. 7. The temperature profile for the harmonic chain, formula (68), for coupling parameter v = 0.05, 
0.2, 1.0 (solid, dotted and dashed lines respectively). 

Eq. (B.5) implies that the value of ^ depends only i — k rather than i + Zc, as before. Physically, 
this symmetry reflects the fact that jj = j independent of the lattice position i in the stationary 
state. In the hmit of large N, Eqs. (B.12,B.13) imply 



2A 



(T+ - T_ 



(70) 



Accordingly, the heat flux is proportional to the temperature difference rather than to the 
gradient as it should be, were the Fourier law to be satisfled. This proves that, as expected, 
homogeneous harmonic chains do not exhibit normal transport properties since the effective 
conductivity k, — jN/ (r+ — TL) oc N, while the bulk conductivity diverges exponentially, since 
the temperature gradient away from the baths is exponentially small. 

For what concerns the dependence of the flux on A, we see that J vanishes both in the limit of 
large and small couplings. The asymptotic expressions 



^ XlksiT^-T^) 



A > a; 
A<a; ' 



(71) 



are thus consistent with the heuristic formula (53) derived in the general case. The maximum 
flux is attained for X/cu = \/^/2, a value that is close to the one observed numerically in the 
nonlinear case (see again Fig. 5). 



Let us conclude this section by recalling that a similar procedure can be adopted to solve the 
problem for heat baths characterized by stochastic elastic collisions. In Ref. [51], it is shown that 
very similar expressions are found also in this case, with only minor quantitative differences in 
the numerical factors. Furthermore, it is worth mentioning the model of self-consistent reservoirs 
introduced in Ref. [53] , that can be solved exactly. 
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Fig. 8. The first, 41st and 100th eigenvector (from top to bottom) in a chain of 130 particles with ran- 
dom masses with an even distribution of 1 and 1/2 values. The increasing localization with increasing 
eigenvalue is transparent. 

4-2 Disordered chains 

While remaining in the realm of harmonic systems, we now consider the role of disorder on 
transport properties. More specifically, we shall consider random-mass (or isotopically disor- 
dered) chains 

mnQn = qn+l " '^Qn + Qn-l ■ (72) 

As we shall see, boundary conditions play a crucial role, but, for the moment, we leave them 
unspecified. Before entering in a more detailed discussion it is worth mentioning the general re- 
sults by Lebowitz and collaborators [27,29]: they showed rigorously that the system approaches 
a unique stationary non-equilibrium state for a large class of heat baths. 

As it is known, the presence of disorder generally induces localization of the normal modes of the 
chain and one may thus expect the latter to behave as a perfect thermal insulator. Nonetheless, 
the actual situation turns out to be much more complicated, depending on boundary conditions 
and on the properties of the thermostats. 

Localization of the eigenmodes - To understand the transport properties it is first useful to recall 
some basic facts about localization. For illustration, let us consider the example of a disordered 
chain with two evenly-distributed types of particles. Some of the numerically computed eigen- 
vectors are shown in Fig. 8. Upon ordering them with increasing eigenfrequencies, a distinct 
difference in their localization properties can be recognized. Indeed, for small frequencies (up- 
per panel of Fig. 8), randomness induces only a relatively weak modulation of the amplitude; 
a partial localization can be recognized in the intermediate panel, while a clear evidence of 
localization is visible only for the high-frequency eigenvector reported in the bottom panel. 
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A rigorous investigation can be performed by applying the transfer matrix approach to the 
eigenvalue equation. After inserting the expression gf„ = v„e''^* in Eq. (72), we obtain 

- mnOu'^Vn = Vn-l ~ 2Vn + Vn+l- (73) 

It is well known that the spectral properties of linear operators involving the discrete Laplacian 
can be determined from a recursive equation for the new variable /?„ = Vn/vn-i- The most 
known example where this approach has been successfully employed is that of Anderson quan- 
tum localization in the tight- binding approximation (see, e.g. [54]) . In the present context, 
Eq. (73) yields 

Rn+i = 2 - m„a;' - (74) 



an equation that can be interpreted as a "discrete time" stochastic equation. The mass rrin 
plays the role of a noise source (with bias), whose strength is gauged by the frequency uj. In 
particular, the inverse localization length 7 is given by 

7=(lni?„), (75) 

while the integrated density of states /(a;) follows from node counting arguments, i.e. = f, 
where / is the fraction of negative Rn values. In Fig. 9 it is shown that / increases hnearly for 
small u and exhibits some irregular fluctuations at larger frequencies. The upper band edge 
(at u ~ 2.8) is easily identifiable as the point above which I{uj) remains constant and equal to 
1. At variance with the standard Anderson problem, where all eigenmodes are exponentially 
localized, here 7 tends to zero for a; — > 0. This can be easily understood from equation (74): 
Since o;^ multiplies the stochastic term, disorder becomes less and less relevant in the small 
frequency limit. In this limit, one can thus resort to a perturbative approach. Let us start 
noticing that for a; = 0, i? = 1 is a marginally stable fixed point of the recursive equation 
(74). For small lu, an intermittent process sets in: after a slow drift driving Rn below one, a 
re-injection to values larger than one occurs and nonlinearity become suddenly relevant. The 
process repeats again and again. By writing R^ — 1 + and expanding in powers of r„, we 
find that the dynamics in the vicinity of it!„ = 1 is described by 

Vn+i ^ rn - r^- (^"^{m) + uj'^5m (76) 

where we have included only the first nonlinear correction and written separately the average 
value of the noise term. In the limit of small cu, this equation can be approximated by the 
Langevin equation 

r = — — uj^{m) -\- cu'^Sm , (77) 
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where, for the sake of simphcity, we have kept the same notations. The corresponding Fokker- 
Planck equation writes 

dP _ d{r'+u;'{m))P ald^P 

dt dr 2 9r2 ^ ^ 



where cr^ = (m^) — (m)^ stands for the variance of the mass distribution. Given the steady 
incoming and outcoming flow, the stationary solution can be obtained by imposing 

{r^ + u;^m))P+^—^C (79) 
z ar 



where C represents the probability flux to be determined by imposing the normalization of the 
probability density P. Notice also that C can be identified with the integrated density of states 
I{uj), since it corresponds to the probability that, at each iterate. P„ is re- injected to the right, 
i.e. the probability of having a node in the eigenvector. In the absence of disorder {am = 0), 



I{uj) = C = 




and, correspondingly. 



1 U\ (m) 



This approximation is already sufficient to reproduce the behavior of I{uj) at small frequencies, 
as in the limit uj the variance of the disorder goes to zero faster than the average value. This 
is confirmed by comparing the dotted line in Fig. 9 (corresponding to the analytic expression 
(80)) with the numerically determined integrated density. 

On the other hand, the above approximation is not accurate enough to determine the local- 
ization length, as disorder is totally disregarded. Indeed, the symmetry of Pq implies that 
7 (r) = 0. By going one step further, we can write P(r) as Pq plus a small perturbation. A 
simple calculation shows that 

uj^a'irJ (m) 

P{r) = Po(r) + , ^2 2 L (82) 



Prom expression (75) for the inverse localization length 7, we find that 



7=(r)-^, for .^0 (83) 
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Fig. 9. Integrated density of states, I{uj) (solid line), and inverse localization length, 7 (dashed line), 
as a function of the frequency u for a chain with mass disorder: the particles have either mass 1 or 
1/2 with equal probability. The dotted line corresponds to the analytic expression (80). In the inset 
the inverse localization length is plotted in doubly logarithmic scales (circles) and compared with the 
theoretical formula (83) - solid line. 
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Fig. 10. Integrated density of states, /(w) (solid line), and inverse localization length, 7 (dashed line) 
for a random-mass chain as in Fig. 9 with the addition of a unit frequency on-site potential. 

an equation derived in Ref. [55] (see Fig. 9). 

If we add a harmonic on-site potential to Eqs. (72), the corresponding scenario becomes analo- 
gous to that of the Id Anderson problem, with all eigen-functions being exponentially localized. 
This is illustrated in Fig. 10, where we have added harmonic springs with unit constant (i.e. a 
force term — g„ acting on the n-th particle) to the chain with random masses considered above. 
The lower band-edge is now strictly bounded away from zero and the inverse localization length 
does not vanish. 

The temperature profile - In order to study the non-equilibrium properties, we need to include 
the coupling with the thermal reservoirs. Here below, we consider Langevin-type heat baths, 
as they allow an analytic treatment, though limited to the weak-coupling regime. The starting 
equation writes 

rrinqn = Qn+i - 2g„ + g„-i + Sni{C+ - Xqi) + SnwiC- - Mn) (84) 
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where for simplicity we have assumed A+ = A_ = A. Although the equations are still linear, 
there is no general method to derive an analytic solution for generic values of the coupling 
constant A. Accordingly, we restrict ourselves to considering the perturbative regime A <C 1. It 
is convenient to introduce the new variable 

Un = Vm^qn , (85) 

which allows rewriting Eq. (84) as 



Un = , 2 h ^ + 

^m„m„+i TTln y/ninmn-l 

— (\/^C+ - Ami) + — (Vm^ (86) 
nil nT'N 

The advantage of this representation is that the operator describing the bulk evolution is 
symmetric and, accordingly, is diagonalized by an orthogonal transformation. In other words, 
upon denoting with the n-th component of the k-th eigenvector, it turns out that J2n ^n^n — 

nj ■ 



Skh and Efc e^e^ = 5., 



With reference to the new variables C/fc = Y,n Un&n ' the equations of motion write as 

k k 

Uk = -colUk - A ^ Ckjilj + + ^^C- (87) 



where —cu^ is the real, negative k-th eigenvalue of the unperturbed evolution operator and 

_e^^e^ . (88) 
mi niN 



Eq. (87) shows that the normal modes are coupled among themselves through the interaction 
with the reservoirs. Standard stochastic calculus applied to the modal energy Ek = (t/| + 
ulU'^)/2 shows that the stationarity condition for the time average {Ek) = implies 

(pk\2 ( k \2 

Ckk{ul) + E CkAUkU,) = + . (89) 



Let us now show that, in the small-coupling limit, this sum is negligible. In fact, from the 
equality 

^ ' = , (90) 

(JjL 
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we find that 

{tlkifh) - ujliUkUh) - Xj2Cjh{tfjUh) = (91) 

j 



By solving this equation together with the symmetric expression obtained by exchanging k and 
h, it is transparent that (UkUh) is proportional to A for A; 7^ h. Accordingly, up to first order in 
A, one has 

m-^(T^^^+TJ^] , (92) 
Ckk \ mi niN J 



As a consequence, by neglecting first order corrections, the local temperature reads 



This is basically the expression derived by Matsuda and Ishii [55] . As a consistency check, one 
can easily verify that if T+ = TL = T, the local temperature Tn is equal to T for all values of n 
(this follows from the normalization condition on the eigenvectors). Furthermore, the profile is 
fiat also when T+ > T_ and the amplitude of all eigenvectors is the same at the two chain-ends: 
in this case, T„ = (T+ + T_)/2. An obvious limiting case is the homogeneous chain. Generally 
speaking, even though the dynamics of a generic disordered chain is statistically invariant under 
left-right symmetry, the same does not hold true for each individual eigenvector. This induces 
some spatial dependence that we shall investigate in the following. 

Visscher [56] challenged equation (93) by arguing that quasi-resonances could be generic enough 
to affect typical realizations of the disorder. In fact, a crucial assumption in the derivation of 
the expression for the temperature profile is that cross-correlations in Eq. (93) are neghgible. 
This is basically correct unless pairs of frequencies are sufficiently close to each other, in which 
case the resonance phenomena should be properly taken into account. Visscher indeed discussed 
particular examples of mass distributions, where a more refined theory is needed. However, as 
long as one is interested in generic realizations, the problem is whether quasi-degeneracies in 
the spectrum are sufficiently frequent to significantly affect the overall vahdity of formula (93). 
In all the cases we have considered this issue turned out to be practically irrelevant. 

Formula (93) does not allow to obtain an analytic form of the profile since it requires the 
knowledge of the eigenvectors and, on the other hand, the locahzation length alone does not 
suffice to predict their amplitude at the boundaries. Therefore numerical diagonalization of 
the Hamiltonian for each different realization of the disorder is required. In Fig. 11, we have 
plotted the stationary temperature profile for a single realization of the disorder versus the 
rescaled lattice position x = n/N. Strong fiuctuations accompany an average decrease of the 
temperature from to TL. 
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Fig. 11. Temperature profile as predicted from Eq. (93) for a given realization of disorder in a chain 
of length = 128, with T+ = 10 and T_ = 5. 




Fig. 12. (a) The disorder-averaged temperature profile as predicted from the formula (93) for differ- 
ent chain lengths (dotted, dashed, dot-dashed and solid curves refer to TV = 64, 128, 256, and 512 
respectively), (b) The variance of the temperature in the same notations as in panel (a). 

This suggests averaging over independent realizations of the disorder in order to better inves- 
tigate the convergence properties with the system size. In Fig. 12a we have plotted the profile 
averaged over 1000 realizations. Upon increasing the chain length, the profile seems to slowly 
attain a linear shape, but sizeable deviations are still present for chains as long as N — 512. 
Such a slow convergence is confirmed in Fig. 12b, where we plotted the sample-to-sample vari- 
ance of the temperature field. Although its asymptotic behavior is even less clear, it is at 
least evident that fluctuations do not vanish in the thermodynamic limit. This is tantamount 
to saying that the temperature profile is not a self-averaging quantity. 

Such difficulties dramatically emerge when performing direct simulations of a disordered chain. 
This issue is of great practical importance also in view of more complex models where analytical 
results are not available. The major problem is represented by the extremely slow convergence 
towards the asymptotic regime that can be explained as follows. Eq. (87) shows that the effective 
coupling of each eigen-mode with the reservoirs is proportional to its square amplitude at the 
extrema. Therefore, all eigenmodes that are localized away from the boundaries can thermalize 
only in astronomically long times. To be more specific, the coupling strength of an eigen-mode 
characterized by an inverse localization length 7 is of order exp(— 7A'"), since it is localized at a 
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Fig. 13. Average temperature profile of random-mass drains for = 16, 32 and 64 (dotted, dashed 
and solid curve, respectively). The coupling constant is A = 1. The average is performed over 1000 
realizations of the disorder each of 5 • 10^ time units. 
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Fig. 14. Behavior of the time average of the temperature in the 5 central sites for A'^ = 16, 32 and 64 
(dotted, dashed and solid curve, respectively). 

distance equal, on the average, to the half of the chain length. This implies that the asymptotic 
profile is attained over times that grow exponentially with N. In other words, the stationary 
state is never reached in the thermodynamic limit. To practically illustrate the issue, we have 
simulated a chain in contact with two stochastic heat baths operating at the same temperature 
T = 5 (see Fig. 13). To further emphasize the slow convergence, all atoms have been initially 
set at rest in their equilibrium positions. Even for relatively short chains (A^ ~ (9(10^) ) almost 
10^ time units do not suffice to fully thermalize the bulk. A more direct way of looking at the 
convergence to the flat temperature proflle is by monitoring the cumulative time average T^, 
performed on the 5 central sites (and averaged also over different realizations of the disorder). 
The data reported in Fig. 4.2, with the choice of a logarithmic scale for the time axis, give an 
idea of the time needed to reach the equilibrium value Tm = 5. 

Heat flux - In the case of stochastic heat baths, like those considered in the previous section, 
one can determine the total heat flux by using Eq. (39). By making use of Eq. (93) we obtain 
[55,56] 

,iX,N)=MT^-T^^ ^^ljf^^ ^^J. (94) 
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where the k-th addendum Jk is naturally interpreted as the contribution of the A;-th mode. As 
intuitively expected, the latter is larger for modes that have larger amplitudes at the bound- 
aries and couple thus more strongly with the reservoirs. This interpretation can be justified 
from Eqs. (87). Indeed, in so far as cross-coupling can be neglected, the dynamics of the k-th 
eigenmode is approximately described by the equation 



k k 

Uk^-u;lUk-XCkkifk + -^^+ + ^^C- . (95) 



Standard stochastic calculus shows that, in the stationary regime, the energy exchanged per 
unit time with the two thermal baths is equal to J^, where coincides with the expression 
imphcitly defined by Eq. (94). 

However, heat transport is characterized by more subtle mechanisms than one could infer from 
this simple picture of independent modes. This is immediately understood if we look at the 
general expression for the local heat fiux, Eq. (27), that, in the case of harmonic chains, reduces 
to 

3 = On) = -{qnqn+i)- (96) 



By expanding and in eigenmodes, this equation can be rewritten as 



On) = - E Wh), (97) 



an expression that, in spite of the explicit presence of the subscript n, is independent of n. The 
interesting point that is made transparent by this formula is that a non-vanishing heat fiux 
is necessarily associated with the existence of correlations among different modes. This is all 
the way more relevant, once we realize that diagonal terms with k — h vanish, being (UkUk) 
the average of the derivative of a bounded function. This observation seems to be in contrast 
with the derivation of Matsuda-Ishii formula itself, that is basically obtained by treating all 
modes as evolving independently of each other. Anyway, we should notice that the heat fiux 
is proportional to A and this is compatible with the existence of weak modal correlations. 
In fact, "velocity- velocity" or "position- velocity" correlations may arise from the fact that all 
eigenmodes are subject to the same noise source (except for a multiplicative factor) and this 
may well induce a sort of synchronization among them. 

The existence of this type of coherence has been numerically investigated and confirmed in 
Ref. [57], where the behavior of a homogeneous harmonic chain has been thoroughly studied. 

Here below, we proceed with our pcrturbative analysis by deriving an analytic expression. Let 
us start by noticing that the equality d{UkUh) / dt = implies that 

{UkUh) = -{UHUk). (98) 
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This antisymmetry property together with the further equahty d{UkUh)/dt — imply that 

AC.. ((^7,^) + {ill)) - H - ^9 = 2A + r_^) (99) 

After replacing the expression of {U^) (see Eq. (92) ) in the above equation, the latter can be 
solved for {UkUh). 

Instead of discussing the general case, we prefer to illustrate the presence of these correlations 
coming back to the simpler case of a homogeneous harmonic chain ( see also Ref. [57] ) . 
In his context, one can, in principle, obtain a a general expression for the correlations by 
transforming Eq. (B.5) (derived for an arbitrary A value) in k space. However, the calculations, 
though straightforward, are rather tedious. Therefore, we limit ourselves to considering the 
weak-coupling limit. The symmetry of the eigenmodes imply that ii S — h — k is an even 
number correlations vanish, while for an odd 5 we have 

rn.^s) - 2 ^^^+"^-^ (100) 

In Fig. 15a we report the numerical results for a chain of length N = 128 with fixed boundary 
conditions and interacting with two thermal baths at temperatures r+ = 75 and TL = 25. 
Apart from the residual statistical fluctuations, a reasonable agreement with expression (100) 
is found upon letting A = 0.056 (approximately equal to the inverse of the average separation 
between consecutive collisions). Fig. 15b shows the results for the same length but a stronger 
coupling strength. The different shape of the curves is a clear indication that higher order terms 
must be taken into account, since the perturbative approach implies that the coupling constant 
acts just as a multiplicative factor. It is anyhow interesting that shape itself is invariant under 
change of 5 as it can be seen in the inset where the three curves are rescaled to their maximum 
value. 

The thermal conductivity - If obtaining an accurate analytic estimate of the heat flux is as 
difficult as determining the temperature profile, we can at least make use of Eq. (94) to deter- 
mine its scaling properties. In fact, since high-frequency eigenmodes are strongly localized, it 
is clear that only the first part of the spectrum contributes significantly to the heat flux. Let 
thus Nf. be the number of modes whose localization length is larger than the sample size N . 
Prom Eqs. (80,83), it follows that 7 ~ a'^I{ujY / {m). Upon writing I — Ne/N and imposing 
^ — 1/N,we find that 

TV, = . (101) 

At this point, it becomes crucial to specify the boundary conditions. Let us first consider the 
case of free ones: the square amplitude of an extended eigenmode in a lattice of size N is of 
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Fig. 15. Modal correlations for a harmonic chain of length A'^ = 128 with fixed boundary conditions. 
Averaging has been performed on a time t = 10^ units. Solid, dashed and dotted curves correspond 
to 5 = 1, 3, and 5, respectively, while the thin lines correspond to the analytic expressions (100). (a) 
refers to a weak-coupling case: the times between consecutive collisions are uniformly distributed in 
the interval [19 -h 21]. (b) corresponds to a strong coupling: collision times distributed in [0.9 ^1.1]. 
The inset contains the same curves, after rescaling to the maximum values. 

the order This implies that the contribution to the heat flux of one of such modes is 

A(T+ — TJ)/N and the heat flux in Eq. (94) can be estimated as 

i^,ee(A,iV) cx A(T+-T_)^-^ . (102) 
As a result, the conductivity diverges as 

Kfree OC \^\fN . (103) 

Cm 



This scaling was first derived in Ref. [55] and later confirmed in Rcf. [58] by means of a different 
approach. On the other hand, for fixed boundary conditions the result is completely different. 
In this case all eigenmodes must vanish for n = and n — N -\- 1. By approximating the 
site-to-site variation of with the wavenumber k/N , we find that the square amplitude of e\ 
and is of order k'^/N^. As a consequence, summing all such addenda up to /c = N^. in Eq. 
(94) yields 

jfU\N) OC A(r+-r_) (^)']^ ■ (104) 

and, accordingly, the thermal conductivity vanishes as 

Kf,^^x(^X^ . (105) 
\crm J VN 

The above estimates give only the leading orders in N. In view of the previously encountered 
strong finite-size effects, it is crucial to check directly the convergence to the asymptotic results. 
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Fig. 16. The effective exponent defined in (106) versus the chain length. The logarithmic derivative has 
been evaluated with finite differences, subsequent points correspond to chain of double length. Circles 

are obtained from the Matsuda-Ishii formula, while diamonds correspond to a direct simulation of a 
chain interacting with stochastic baths operating at T+ = 10 and T- = 5, respectively. The collision 
times were uniformly distributed in the range [1-^2]. 

To this aim it is convenient to compute the effective exponent 



The results are shown in Fig. 16 for the case of fixed boundaries. For weak coupling, the 
conductivity has been evaluated by numerically computing the eigenvectors and averaging the 
Matsuda-Ishii formula (94) over 1000 realizations of the disorder. The asymptotic regime a = 
— 1/2 is approached very slowly (see the circles): one should consider N values much greater 
than 10^. Similar results are found at stronger coupling by directly simulating chains that 
interact with stochastic baths. The data (diamonds in Fig. 16) suggest that a relatively strong 
coupling reduces the amplitude of finite-size corrections. Finally, it is important to realize that 
the small coupling of the exponentially localized modes with the thermal baths docs not cause 
any problem to the temporal convergence of j(A, A^), since independently of whether such modes 
have reached their stationary state, their contribution to the total heat flux is anyhow negligible. 

In summary, not only boundary conditions affect the scaling behavior of k, but they give rise 
to qualitatively different scenarios: for free boundaries, disordered harmonic chains exhibit an 
anomalous conductivity as it diverges in the thermodynamic limit. On the contrary, a disordered 
chain with fixed boundaries behaves as good insulator! This latter scenario is brought to an 
extremum if we add an on-site potential. In fact, we have already mentioned that all eigen- 
functions become exponentially localized and this implies that conductivity is exponentially 
small in N. This is again very much reminiscent of the electrical conductivity of the Anderson 
problem. 

Dhar [60] went even further and showed how the scaling behavior of the conductivity with the 
system size depends also on the spectral properties of the heat baths. More precisely, if k oc A^*^ 
then the exponent a is determined by the low-frequency behavior of the noise spectrum. This 
implies that a suitable choice of the latter can even lead to a finite conductivity! Such a scenario 
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Fig. 17. The heat flux (scaled to the maximum value) versus the wavenumber k (scaled to the position 
of the maximum) for different chain-lengths (dotted, dashed, dot-dashed and solid curves refer to 
N = 32, 64, 128 and 256, respectively. The same quantities are reported in doubly logarithmic scales 
in the inset, to show the quadratic growth for small k. 



is less unphysical than it may appear at a first glance. Integrability of the motion implies that 
the only scattering mechanism that determines the heat resistance is the interaction with the 
baths. It is therefore reasonable that the actual way in which the latter transfer energy among 
the modes plays a crucial role. 

Modal fluxes - We conclude this section with a discussion of the modal heat fluxes Jk as 
defined from Matsuda-Ishii formula (94). Besides providing a finer verification of the latter, the 
analysis is useful in understanding the individual contributions of each "channel" to the heat 
transport. The spectra Jk obtained for different chain lengths are reported in Fig. 17. They 
have been scaled each to the maximum Jm, while k has been scaled to the wavenumber /cm of 
the maximum itself. In practice, this is asymptotically equivalent to scaling the vertical axis 
by a factor N'^ and the horizontal axis by 1/\/N. We have preferred to adopt this strategy in 
order to possibly get rid of the strong finite-size corrections revealed by the previous analysis. 
Indeed, the good data collapse (except for the right tail) noticeable in Fig. 17 is very suggestive 
of the existence of an asymptotic spectrum. Additionally, the quadratic growth predicted for 
fixed b.c. is much more clear (see the inset) than one could have expected from the scaling 
behavior of the total heat flux for the same chain lengths. 

Besides looking at the average heat fluxes Jk, we have studied their sample-to-sample fluctua- 
tions, by computing the variance aj. The relative variance plotted versus the scaled (as in the 
previous figure) wavenumber indicate that fluctuations are independent of the chain length (see 
the almost overlapping curves in Fig. 18): this means that Jk is not a self-averaging quantity 
in the thermodynamic limit and this holds true even for low /c-wave-numbers, whose behavior 
results, in principle, from a spatial averaging over increasingly long spatial scales. More pre- 
cisely, we observe that the variance increases exponentially with k, starting approximately from 
0.26 for the longest wavelength. Moreover, we have compared the theoretical results with non- 
equilibrium simulations with stochastic heat baths and collision times distributed uniformly in 
the interval [30 -j- 60] . The data plotted in Fig. 19 reveal a very good agreement over various 
orders of magnitude if the coupling constant is set equal to 1/54. 
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Fig. 18. The relative standard deviation versus the rescaled wavenumber k reveals a clear exponential 
growth. 



Fig. 19. Comparison between numerical results and the theoretical formula for the heat flux (solid 
curve) (see eq.(94) ). The agreement is fairly good for both N = 32 (squares) and A'^ = 64 (diamonds). 

5 Linear response theory 

5.1 The Boltzmann-Peierls equation 

In this section, we briefly review the "traditional" approaches to the determination of thermal 
conductivity. Although they have a major importance in solid state applications, here, we limit 
ourselves to a very sketchy discussion, since it is sufficient to point out only those general issues 
that are of interest for our purposes. In this respect, our presentation is inspired by the review 
article of Jackson [61]. 

The most elementary picture of heat conductivity is based on the analogy with kinetic theory 

of gases where k, = Cvgi/S, C being the heat capacity, Vs the sound velocity and C, the mean free 
path. In a lattice, one can imagine to replace particles with normal modes, but it is, of course, 
necessary to take into account that the latter have different group velocities, fk = duj/dk, 
depending on their wavenumber. Accordingly, the above expression for k generalizes to 
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where we have introduced the relaxation time Tk = ^k/'t'k that can be determined by phcnomcno- 
logically inchiding all possible scattering mechanisms (anharmonicity, impurities, boundary ef- 
fects, electrons etc.) that must be computed in some independent way. 

A less heuristic derivation of the above formula is obtained by solving the Boltzmann equation 
in the relaxation time approximation [62]. In 1929 R. Peierls proposed his celebrated theoretical 
approach based on the Boltzmann equation. The main idea is again taken from kinetic the- 
ory: lattice vibrations responsible for heat transport can be described as an interacting gas of 
phonons [2]. Accordingly, one can introduce the time-dependent distribution function Nk{x,t) 
of phonons with wavenumber /c in a macroscopically small volume around x. ^ If we further 
limit ourselves to considering only the cubic term in the interaction potential (three-phonon 
processes), the kinetic equations are of the form 

dNk , ON, 



dt 



+ ^k-^ = J J dk'dk"{[NkNk>Nk>'-{Nk + l){Nk> + l)Nk>>]Wkk'k" 

+ ^[Nk{Nk> + l){Nk'> + l)-{Nk + l)Nk^N„.]Wkk'k"} , (108) 



where the transition probability Wkk'k" is basically obtained from the Fermi's golden rule. The 
r.h.s. is the collision integral, i.e. the difference between the number of processes (per unit 
time) that either increase or decrease the number of phonons in the state k. These nonlinear 
integro-differential equations are clearly impossible to solve in general. An approximate solution 
is obtained in the limit of small applied gradients, i.e. by looking for small perturbations of the 
equilibrium distribution — N^^ + 5Nk, where — {exp{hwk/kBT) — 1)~^. This allows to 
write a linearized kinetic equation [2,48] that, in the stationary case, is of the form 



where I is the linearized coUision integral, which is a linear functional of the SN^s. The Tk are 
thus determined as the eigenvalues of the problem. 

Anyway, some useful information about thermal conductivity can be obtained by looking di- 
rectly at the dynamics in Fourier space. Whenever third and fourth order terms are present in 
the equations of motion (as in in the FPU model (4)), one can write 

Qk = -^iQk - X! ^kkik2QkiQk2 - X) ^kklkikaQkiQkiQka ■ (HO) 
ki,k2 k\,k2,k3 

Accordingly, the harmonic part of the flux given by Eq. (33) satisfies the dynamical equation 

Jh =^-YY^ [-^ki^k + Vk^UJki+Vk2l^k2]y-kkik2QkQkiQk2 
kkik2 

^ Here, for simplicity, we are referring to a one-dimensional ordered crystal. 
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+ [-'"ki^k + Vkii^ki + Vk2i^k2 + Vka^^ka] V-kkrk2k3QkQkiQk2Qk3- (HI) 

—kkik2k3 

As expected, this implies that Jh is a constant of motion in the harmonic case. In a perfect 
lattice, one has the selection rules (remember that the mode indices range between —N/2 + 1 
and 



^-feUv^O for -k + h + k2^0,±N 

^itWsT^O for -k + h + k, + h^O,±N. (112) 

Peierls observed that if there is no dispersion, i.e. ujk = Vs\k\, both three- and four-phonon 
contributions to (HI) vanish when the sums in (112) are equal to zero. Therefore, the only 
finite contributions to Jh are those arising from the so-called Umklapp processes corresponding 
to the above sums being equal to ±A^. 

Besides the above general considerations, there are some specific comments regarding the role 
of dimensionality that can be drawn in the framework of perturbative theories. Indeed, by 
evaluating the r.h.s of (111) to lowest order (i.e. by replacing Qk with the harmonic solution 
(35)) and averaging out the fast oscillations, one is left with the leading resonant terms that 
satisfy additional conditions like 

- OUk + OUki + 0Uk2 ^ - OUk + OUki + OUk^ + OUka ^ (113) 



etc. Thus there is a big difference between three and four phonon processes in 1-dimension, as 
in the former case the first of conditions (112) and (113) cannot be simultaneously satisfied 
(see Ref . [63] for some numerical results) . The net result of this argument due to Peierls is that 
the lowest-order contributions to thermal resistance in 1-dimcnsion arises from four phonon 
Umklapp processes. Of course, in higher dimensions, the situation is different, because the 
energy and momentum constraints can be satisfied also by three-phonon processes due to the 
existence of different (longitudinal and transverse) branches of the frequency spectrum. Let us 
finally notice that, although - in the spirit of a perturbative calculation - only the harmonic 
component of the flux Jfj has been considered, no basic difference arises upon including also 
the nonlinear component [61]. 

Similar conclusions can be drawn from the analysis of the high-temperature limit of the Boltz- 
mann equation, following a standard argument originally due to Pomeranchuk (see for exam- 
ple Chapter VII in Ref. [48]). In this limit, with reference to processes involving three long- 
wavelength phonons of wavenumber k, k', k", the transition probability Wkk'k" ~ kk'k" ~ A;^, 
whereby N*^ ^ ksT/hiOk ~ 1/k. From this consideration and from Eq. (108), one can esti- 
mate the linearized collision integral appearing in Eq. (109) as X ~ 6Nk/Tk ~ k'^^^dNk, in 
dimension d. Accordingly, the solution of Eq. (109) diverges as SNk ~ A;"^^''"'^). Finally, using 
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the expression for the heat flux (33) with {Ek) — huJkNk yields 

^^^^^ 

This means that the contribution of such processes would lead to a thermal conductivity diverg- 
ing like in any dimension. In order to avoid this divergence it is therefore necessary that long 
wavelength phonons are scattered by short- wavelength ones. Following Ref . [48] , let us consider 
the process in which a short-wavelength phonon of index k annihilates into two phonons of 
index k' and k — k' — N with k <^ N and k <^ k' . The condition (113) requires ujk — ^^k-k' +<^A;' 
that, in turn can be satisfied only for Vk — Vs- More generally, one can show that, in the absence 
of degeneration points in the spectrum, the condition to be fulfilled is that the group velocity 
of short wavelength phonons is larger than the sound velocity i.e. \vk\ > Vg- Once again this 
constraint cannot be satisfied in onc-dimensional homogeneous chains and one would conclude 
that a finite conductivity can be possibly established only by means of higher-order processes. 

The Boltzmann-Peierls approach is certainly one of the milestones of the theory of thermal 
transport in solids. Nonetheless, it is important to recall that its derivation is essentially based 
on second-order perturbation theory (through the collision kernel Wkk'k", which is evaluated 
by means of Fermi's golden rule) and involves the use of random phase approximation among 
the phonon modes, which is certainly less appealing than the Stosszahlansatz originally intro- 
duced by Boltzmann for molecular collisions. It is however remarkable to notice how classical 
perturbative approaches are able to predict some peculiarities in low-dimensional anharmonic 
lattices. 

5.2 The Green-Kubo formula 

The other major tool, commonly used when dealing with transport processes, is linear response 
theory. At variance with the response to mechanical perturbations (e.g. an external electric 
field), heat conduction is a process driven by boundary forces. Therefore, a conceptual difficulty 
arises, since there is no explicit small term in the Hamiltonian to be used as an expansion 
parameter. This difficulty can be overcome at the price of a stronger assumption, namely that 
local equilibrium holds. The hypothesis looks physically reasonable, but it is far from being 
rigorously based even in simple mathematical models and it has been often devised as one of 
the weak points in the foundation of the whole theory. If local equilibrium holds, a temperature 
field T(x) can be defined accordingly, thus allowing to introduce a non-equilibrium distribution 
function 

p = Z-^exp J dx/5(x)/i(x)^ , (115) 

where /i(x) is the Hamiltonian density, while Z is the partition function. By now assuming that 
the deviations from global equilibrium are small, we can write /3(x) = /9(1 — AT'(x)/T) and 
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thus 

p = exp + , (116) 

where Ti' is the perturbative Hamiltonian 

n' j cixAT(x)/i(x). (117) 

It is therefore possible to proceed with a perturbative expansion, obtaining the well known 
Green-Kubo formula that in the classical case reads [25] 

1 \ 

^GK = TT^ / dr hm V-i(J(r)J(0)) (118) 



where J is the total heat flux defined by Eq. (31) and kgk should, more properly, be a tensor. 
However, in the simple case of isotropic homogeneous solids made of atoms placed on a reg- 
ular hyper-cubic lattice, the thermal conductivity tensor has a diagonal representation, with 
equal non-zero components: upon these assumptions, in dimension d, kqk reduces to the scalar 
quantity 

As often stated, Eq. (118) relates the non-equilibrium transport coefficient to the fluctuations of 
a system at equilibrium. It has to be reminded that its rigorous mathematical foundation is still 
lacking [3]. Besides this, there arc several subtleties connected with a correct implementation 
of this formula. First of all, one should notice that the infinite-volume limit should be taken 
before the long-time limit, in order to avoid the problem of Poincare recurrences. This is a 
particularly delicate matter whenever a slow decay of correlations is present. 

The next issue concerns the meaning of the ensemble average (•). In the derivation d la Kubo, 
it denotes a canonical average, while the formally identical expression obtained by Green refers 
to the micro-canonical ensemble. In this latter case, if the total momentum P is conserved, it 
has to be set equal to zero, otherwise (J) 7^ and the integral in Eq. (118) would trivially 
diverge. Alternatively, as observed in [3] , one may compute the truncated correlation functions 
(J(i)J(0))T = (J(i)J(O)) - (J)^ for any Py^O.^ 

Another way to see the problem was pointed out by Green himself [65]. He noticed that the 
microscopic expression of the heat flux to be employed in Eq. (118), depends on the chosen 
ensemble. Indeed, he showed that while Eq. (31) is the correct expression in the micro-canonical 

^ Overlooking this point may lead to some confusion as in Ref. [64]. 
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case, a "counter-term" must be subtracted in the grand-canonical case. This is readily seen by 
letting Xj — > Xj — V where v is the velocity of the center of mass. Up to terms of order less than 



where E is the average energy, p the pressure and H the total enthalpy of the system. The 
micro-canonical and grand-canonical ensembles give the same results provided tliat the micro- 
canonical energy density is chosen to correspond to the canonical temperature T. The reason for 
the different expressions to be used is that the same observable has different time-correlations 
in the different ensembles. 

5.3 Mode-coupling theory 

Despite the conceptual difficulties lying behind the derivation of Eq. (118), this formula provides 
a well defined prescription for determining the thermal conductivity kgk from the current- 
current correlation function at equilibrium. We claim that an effective method for estimating 
this correlation function is provided by the well known mode-coupling theory (MCT), intro- 
duced some decades ago to approach the problem of long-time tails in fluids [5] . Since a rigorous 
proof of this statement is still lacking, we prefer to illustrate first some simple arguments to 
support the claim. Afterwards, we introduce MCT in the simple context of the FPU model with 
cubic nonlinearity. Finally we briefly recall the major quantitative results in various dimensions. 

According to the classical perturbative approach outlined in the flrst section of this chapter, 
the time scale for the relaxation process towards the stationary state can be determined by 
linearizing the collision operator in Eq. (108). However, this may be insufficient, because of 
the possible existence of subtle dynamical correlations that escape the predicting ability of a 
perturbative approach: see, for instance, the necessity to invoke higher order processes at high 
temperatures, or the well known existence of slow relaxation processes at low temperatures 
when the dynamics is almost integrable. 

A more powerful approach can be built starting from the observation that in solids, like in 
fiuids, the slowest processes arise from the "diffusion" of conserved quantities (such as energy 
and momentum). Actually, macroscopic conservation laws necessarily imply the existence of a 
hydrodynamic behavior, dominated by the time scales l/'~fk associated with the dynamics of 
long-wavelength, i.e. low-k, modes. It is then crucial to observe that the damping factor 7^ is 
expected to vanish in the limit /c — > both in crystals characterized by the existence of an 
acoustic band, as well as in fiuids (e.g., hard spheres interacting via short-range potentials). In 
fact, as k = modes correspond to exactly conserved quantities, long-wavelength ones must, 
by continuity reasons, be characterized by a slow dynamics. It is then of primary interest to 
notice that this is true independently of the strength of the perturbative terms. Accordingly, 
7fe may be very small also for strong nonlinearities, when no standard perturbative approach 
can be meaningfully implemented. As a matter of fact, numerical studies of several models of 





{E + pV)v = Jmic-Hv 



(120) 
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anharmonic crystals with confining nearest-neighbor interactions (hke the FPU model) nicely 
confirm this scenario [66] also at high temperatures, when high-/c modes behave like "thermal" 
variables, rapidly relaxing to equilibrium. ^ 

According to the previous discussion, the long-time behavior of any current-current correlation 
function depends on how 7^ — in the limit A; — > 0. In particular, if the damping factor 
vanishes too rapidly with /c — 0, the temporal decay of the heat-fiux correlation function may 
be so slow that the integral in Eq. (118) diverges. In general this effect ought to depend also 
on the space dimension. Indeed, almost conserved modes propagating with the sound velocity 
through the lattice are expected to propagate more efficiently in low than in high dimensions, 
where the presence of transverse modes favors collision mechanisms. Actually, a well defined, 
i.e. finite, transport coefficient in anharmonic solids should emerge from an efficient dissipation 
of the energy of sound waves. In this sense, it is worth recalling that Fourier law follows from 
the assumption that the temperature field obeys a difTusive equation. 

As already anticipated, we illustrate an apphcation of MCT to the FPU model with cubic 
nonlinearities. The procedure is an extension of linear response theory and represents a first 
step towards the construction of a formal approach for the description of transport properties 
in models of Id solids. Moreover, it provides the theoretical background for arguing that the 
same features should be observed in all models where nonlinear effects can be ascribed to the 
two leading algebraic terms. 

In the framework of linear response theory, the dynamics of slow modes is described by gener- 
alized Langevin equations. These are linear stochastic equations with memory terms and are 
usually derived with the projection method introduced independently by Mori and Zwanzig 

[25]. To illustrate this in the present context, let us consider a one-dimensional chain like (2) 
with periodic boundary conditions. The equations of motion for the normal coordinates (32) 
can be written as (see also Eq. (110)) 

Qk = -ojlQk + J'k (121) 



where J^k accounts for mutual interactions among the modes while uj^ denotes the normal-mode 
frequency. One can then define a projection operator P, acting on the generic scalar observable 
O as 



k 



(OQl) 



Qk + 



(OQl) 



Qk 



(122) 



' In the high-energy regime, the same models are known to exhibit a strongly chaotic and convincingly 
ergodic behavior. The time-scale separation between low-A; and high-fc modes seems to contradict this 
statement. This is not the case, since it has to be noticed that low-k modes, although playing a major 
role in transport phenomena, are a negligible fraction of the spectrum in the thermodynamic limit. 
Accordingly, equilibrium properties are dominated by thermal modes and "hydrodynamic" deviations 
from ergodicity can be detected only as higher order corrections. 
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Due to the conservation law of total momentum, we expect that the slow dynamics should be 
associated with the long-wavelength Fourier modes Qk with \k\ ^ N/2. Moreover, translational 
invariance implies that each mode is uncorrelated from the others, so that we can consider each 
mode separately. The corresponding projected equations of motion (that are still exact) read 
as [68] 

t 

Qk + jTk{t-s)Qk{s)ds + CjlQk = Rk , (123) 



where the random force Rk{t) — {l—V)Qk is related to the memory kernel by the fluctuation- 
dissipation theorem 

rk{t)^P{Rkmi{0)). (124) 



The first effect of nonlinearities is to induce a temperature-dependent renormalization of the 
dispersion relation 



This amounts to renormalizing the sound velocity from the "bare" value Vg to Vs — VgV^ + Oi. ^ 
A straightforward consequence of Eq. (123) is that the normalized correlation function Gkit) — 
(5Cjl{Qk{t)Q%{Q)) (^fe(O) — 1) obeys the equation of motion 

t 

Gkit) + J Tkit - s)gk{s) ds + ulSkit) = . (126) 





Up to here we performed an exact but formal manipulation of the equations of motion. The 
crucial point is the explicit computation of the memory kernel rfe(i). MCT is an approximate, 
self-consistent method for obtaining such an expression in terms of Gkit)- A first conceptual 
difficulty of the projection approach is that R^ does not evolve with the full Liouvillean operator. 
One can bypass the problem with the replacement [5] 

{Rkmim - (Mtmm, (127) 

whose validity is based on the implicit hypothesis that the slow terms possibly contained in Jvc(^) 
can be neglected in the thermodynamic limit. A second simplification amounts to factorizing 
multiple correlations. For example, in the case of a quadratic force, one obtains [68] 

^k{t) « ^J2^>''it)Sk-k'{t) ■ (128) 

k' 

^ Notice that for T — (/3 ^ oo), a{f5) — >■ as the integrals in Eq. (125) reduce to Gaussian integrals. 
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This approximate expression of the memory kernel constitTites, together with Eq. (126), a closed 
system of equations for that has to be solved self-consistently by introducing the Laplace 
transform of F^, 

oo 

rk{z) = J e-'''rk{t)dt, (129) 



and, analogously, Qk{z). One finds that (with ^^(0) = 0) 

iz + Tkjz) 

As long as dissipation is small enough, Qj^ has the form 

gk{t) ~ exp(iAfct), (131) 

where the pole of the above transform is approximately given by 

A. = ±^. + i^^, (132) 



which can be regarded as a generalized dispersion relation. The imaginary part 7^ of repre- 
sents the effective relaxation rate of each Fourier mode as a consequence of its interaction with 
all the other modes. 

As discussed in Ref. [5], in one dimension, the self- consistent calculation predicts a singularity of 
the memory function at z = 0: T{z, q) ~ z'^^^q^. Substituting this result into the approximate 
dispersion relation (132) yields the non-analytic dependence of the relaxation rate in the limit 
of small wave-numbers [67] 

7(g) « q'/'. (133) 

The scaling behavior (133) has been confirmed in molecular-dynamics simulations performed 
at equilibrium for the FPU model [68] in a significative range of chain sizes, (see Fig. 20). The 
above discussion can be generahzed to the case in which a cubic force is present [68] . Although 
the dependence of the relaxation rate on the temperature depends on the specific form of the 
anharmonic potential, the rate q is expected to be the same for all one- dimensional models, 
where the theory applies. Application to 2d lattices has been also discussed [69]. 

More in general (see again [67]), the following dependences are found for d = 1 3, 
X{q) c^cq- idq^'^ + ... {d^l) 
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Y(q) - 




Fig. 20. The wavenumber dependence of the relaxation rates 7(g) for the quartic FPU potential (4) 
54£ = 8.8 . All the points were obtained from the initial decay of the envelope of Qi for increasing 
values of AT up to iV = 2048. The dashed line is a power-law fit q^-^"^. 



A(g) ~ cgr-icV In gf + ... (d = 2) 
X{q) ^cq- ic'q^ + ... (d = 3) 



(134) 



where c' is a suitable multiplicative factor. Direct numerical evidence of the peculiar behavior 
of one-dimensional systems has been found also in fluids [70,71]. 

In order to understand the consequences of the above reasoning for transport phenomena, it is 
convenient to look at the dynamics of Fourier modes. If we assume that memory effects are all 
contained in a mode-dependent relaxation-time, Eq. (123) effectively reduces to its Markovian 
limit, 

Qk + ikQk + ojlQk = Rk (135) 



(referred to the finite- length case), where the random force is well approximated by a Gaussian 
white process 

{R,{t)Rl{t')) = ^6{t-t') , (136) 



and, for the sake of simplicity, we have neglected the small frequency shift possibly arising from 
the solution of (132). At this level of approximation, the physical implications of the mode- 
mode interactions are contained in the dispersion relation 7^. Eq. (135) provides a convincing 
description of the numerical results reported in Refs. [66,68], where the Fourier-mode dynamics 
was studied for the FPU-/3 model. 
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Let us now split V into its harmonic and anharmonic parts and consequently write the flux 
(30) as J = Jh + Ja- For a strongly anharmonic system, like the one we consider here, we do 
not expect Ja to be negligible. Nevertheless, in the spirit of Section 5.1, one can argue that the 
two terms exhibit the same leading asymptotic behavior, so that we can restrict ourselves to 
considering the autocorrelation of Jh alone. This hypothesis has been also successfully tested 
in simulations [72]. 

It is convenient to rewrite expression (36) as 

Jh = Y.Vk{Ek-{Ek)) = Y.Vk5Ek (137) 

k k 



(compare also with Eq. (114) where 6Ek are the energy fluctuations of each mode). Notice also 
that, in view of the renormalization of the frequencies, E^ should now be defined on the basis 
of the transformations (35) with cuk replaced by cuk- Prom Eq. (135), one expects that, for small 
7fe, energy fluctuations satisfy the Langevin equation 

5Ek = -jkSEk + R'k , (138) 



with no oscillation for 6Ek. In such a limit, i?^ is well approximated by a Gaussian and delta- 
correlated random process and {{SEkY) = k'^T'^- For large N, we obtain [72] 



n/a 

Na 



{Mt)JHm oc Y.^,m6E,f)e-^>'' = -^klT' J dqv'{q)e-^^^^' . (139) 



— 7r/a 



Since the integral is dominated by the low-^ contribution at large times, we can estimate the 
long-time behavior of Eq.(139) by letting c{q) ~ Vg and extending the integration to inflnity. 
Furthermore, in accordance with Eq. (133), we let ^{q) — c'q^, obtaining 



{Mt)jHm oc 



)lklT^Na 



1 + o (t-l) 



(140) 



This result can be generalized to derive the following long-time behavior for the heat-flux 
correlation function. 



(j(t)j(o)) - r'/^ {d = 1) 

(J(t)j(o)) (d = 2) (141) 

(J(t) J(0)) ~ r=^/' (d = 3). 

The knowledge of the asymptotic behavior of (J(t) J(0)) now allows determining the dependence 
of K on N. In fact, upon restricting the integral in Eq. (118) to times smaller than the typical 
transit time Na/vg, one obtains. 
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id -I) 
{d = 2) 
(d = 3) 



K ~ finite 



(142) 



It is worth stressing again that these results can be derived without making any exphcit reference 
to the details of the interactions among atoms in the lattice. In practice, this is tantamount to 
stating that all models characterized by short range interactions and momentum conservation 
should exhibit the same kind of anomalous behavior (for < 3). In the next chapter we shall 
see that this is not completely correct, as for non confining potentials, a normal conductivity 
is found already in one dimension. In all other cases, in spite of the intrinsic approximations 
contained in the MCT, the predicted scaling behavior of k agrees with the numerical estimates 
obtained from both equilibrium and non-equilibrium simulations. This is not surprising if one 
considers that, at relatively high energies, the time scale associated with high-fc modes is well 
separated from the hydrodynamic ones, so that the leading term predicted by MCT is likely to 
contain all the relevant information already in relatively small systems. 



6 Anharmonic chains with momentum-conserving potentials 

6.1 Early results 

This section is devoted to a historical review of molecular-dynamics studies of thermal con- 
duction in the class of models (2). The first simulations date back to the pioneering work of 
Payton, Rich and Visscher [32] and to the contribution of Jackson, Pasta and Waters [73]. In 
both cases, the Authors performed non-equilibrium studies of the FPU model (4) with coupling 
constants g2, g^, and chosen in such a way to represent the leading terms of the expansion of 
the Lennard- Jones potential (3). In order to study the effect of impurities in the crystal, either 
a disordered binary mixture of masses [32] or random nonlinear coupling constants [73] were 
considered. Ironically enough, those very first computer studies attacked the problem from the 
most difficult side. In fact, even before the effect of disorder was fully understood in harmonic 
chains, they studied systems where anharmonicity and disorder are simultaneously present. 
Nevertheless, those early works have at least the merit to have showed how the interplay of 
the two ingredients can lead to unexpected results that, in our opinion, are still far from being 
fully understood. Indeed, Ref.[32] revealed that the simple perturbative picture in which anhar- 
monicity and impurities provide two independent (and thus additive) scattering mechanisms 
does not hold. More precisely, the Authors found even cases in which anharmonicity enhances 
thermal conductivity. A qualitative explanation was put forward by claiming that anharmonic 
coupling induces an energy exchange between the localized modes, thus leading to an increase 
of the heat flux. 

The limited computer resources available at that time prevented, however, addressing the issue 
whether the combined effect of anharmonicity and disorder can lead to a finite conductivity. On 
the other hand, it was payed attention at the temperature profile T{x), noticing irregularities 
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that depended on the reahzation of the disorder. While we know that T{x) is not a self- averaging 
observable of disordered harmonic chains, it has not yet been clarified whether the dependence 
on the realization of the disorder persists over long enough time scales in anharmonic chains as 
well. 

Additional questions that have been investigated concern the concentration of impurities. Be- 
sides the obvious finding that disorder reduces the value of heat conductivity (for fixed finite- 
chain length), it was noticed an asymmetric behavior between the case of a few heavy atoms 
randomly added to an otherwise homogeneous light-atom chain and its converse. The smaller 
values of the conductivity observed in the former cases were traced back to the larger number 
of locahzed modes [32] . 

Having recognized the difficulty of simultaneously coping with the effects of nonlinearity and 
disorder, we now turn to the simpler case of anharmonic homogeneous chains. Some early 
work in this direction was performed by Nakazawa [74] who considered equal-masses FPU 
and Lennard- Jones chains composed of 30 particles and coupled with Langevin baths at their 
boundaries. This setup required the integration of a set of stochastic differential equations, a 
task that was admittedly unfeasible with the computer resources available at that time. As a 
consequence, several attempts of designing artificial but easy-to-simulate models followed these 
first studies. Some examples are reviewed in Ref. [75] ^ Let us mention among them the case of 
the harmonic hard-rod potential, i.e. a harmonic well delimited by an infinite barrier located 
at a given distance from the equilibrium position. A diverging conductivity was observed with 
a method akin to Green-Kubo one [75] . 

The long period of time (almost a decade) during which the problem was practically forgotten 
signals perhaps the frustration encountered in the search for the minimal and general require- 
ments for building simple Id models with good transport properties. In the mid eighties several 
authors got again interested in the problem, being able to perform non-equilibrium simulations 
of chains with smooth inter-particle potentials and a few hundreds of particles. In particular, 
a good deal of work was devoted to reconsidering the FPU model[76,77] and to studying the 
diatomic Toda chain [78,79,80,81] whose Hamiltonian is of the type (2), with 



and mi being a sequence of alternating light and heavy masses with given ratio r. At variance 
with the homogeneous case r = 1, this model is no longer integrable at equilibrium and can 
be thus considered as a meaningful candidate for testing the vahdity of the Fourier law. Notice 
that at variance with homogeneous models, (143) admits also an optical branch in the harmonic 
limit. 

Out of the many conflicting results, Mareschal and Amellal [82] recognized that fluctuations of 
the heat current display peculiar features that are normally absent in higher dimensions. More 

^ Those attempts eventually led to the invention of the so-called ding-a-ling model described in the 
following Section as it belongs to the different class of chains with external substrates. 
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Fig. 21. Temperature profiles for the FPU/? model for r+ = 10 and T_ = 8. Panels (a) and (b) refer 
to stochastic reservoirs (acting through the randomization of the velocity at random times uniformly 
distributed in the interval [1,2]). Dotted, dashed, dot-dashed, and solid lines correspond to = 128, 
256, 512, and 1024, respectively. Panels (c) and (d) refer to Nose-Hoover thermostats with = 1. In 
this case, dotted, dashed, dot-dashed, and solid lines correspond to iV = 32, 64, 128, and 256. In (a) 
and (c), fixed b.c. are imposed, while free b.c. are imposed in (b) and (d). 

precisely, they computed the equilibrium autocorrelation function of the heat flux (the Green- 
Kubo integrand) for a Lennard- Jones chain of 200 particles. On a qualitative level, they noticed 
that the initial fast decay was followed by a very slow convergence to zero. This is consistent with 
the possibility that such a long-time tail be responsible for a diverging transport coefficient, 
in close analogy with what happens in low-dimensional fluids [5]. Additionally, in Ref. [82] 
it was checked the robustness of this feature against the introduction of further, extrinsic, 
scattering mechanisms. Indeed, the time tails survive the addition of a moderate fraction of 
impurities, introduced as either mass defects or variable interaction potentials. For instance, it 
was considered the case where every fourth particle has a purely repulsive interaction with its 
neighbors, of the type 



A slow decay was finally observed also in the presence of an external sinusoidal field (akin to 
the Prenkel-Kontorova substrate potential considered in the following). It has, however, to be 
recognized that this last observation contrasts the recent results obtained for several models 
with on-site forces (see Sec. 7). 

6.2 Divergence of heat conductivity 

It took almost another decade for the first systematic studies on the size dependence of the 
conductivity to appear. An extensive series of non-equilibrium simulations were performed for 
the FPU chain with quadratic and quartic [83,84,72] or cubic [63] interaction potential as well 
as for the diatomic Toda one [85,86]. 

Some evidence of anomalous transport properties is already given by the temperature profiles. 
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Model 


Reference 


a (NEMD) 


a (GK) 


FPU-/? 


[84,72] 


0.37 


0.37 


FPU-a 


[63] 


<0.44 


- 


Diatomic FPU r=2 


[86] 


0.43 


compatible 


Diatomic Toda r=2 


[85] 


0.35-0.37 


0.35 




[86] 


0.39 


compatible 


Diatomic Toda r=8 


[86] 


0.44 


compatible 


Diatomic hard points 


[85] 


0.35 





Table 1 

The estimated exponent a of divergence of the conductivity with size N , as obtained from both non- 
equilibrium molecular dynamics (NEMD) simulations and through Green-Kubo (GK) equilibrium 
studies. Only the significative digits are reported as given in the quoted References. 

While a fairly linear shape is obtained for free boundary conditions (see Fig. 21b,d), as pre- 
dicted from the Fourier law, strong deviations are observed for fixed boundary conditions (see 
Fig. 21a,c). More important, such deviations persist upon increasing the chain length: the nice 
overlap observed in panel (c) indicates that the asymptotic temperature gradient is definitely 
non uniform. Altogether, this scenario is suggestive of the existence of long-range effects. 

Before discussing the divergence of let us remark that, as anticipated in Chap. 3, a much 
smaller boundary-resistance and thus smaller finite-size corrections are found in Nose-Hoover 
thermostats than in stochastic ones (see Fig. 21). 

As a result of the recent numerical studies, one can now safely claim that the conductivity of 
long but finite chains diverges as 

k{N) oc AT" (145) 

In Table I we compare the available estimates of the exponent a determined by different authors 
in various models. The numerical values range between 0.35 and 0.44, suggesting a nontrivial 
universal behavior. It is also remarkable to notice the overall consistency among the results 
obtained with different thermostat schemes (ranging from deterministic to stochastic ones). 

In order to better appreciate the quality of the divergence rate that can be numerically obtained, 
in Fig. 22 we have plotted the finite-length conductivity k{N) = JN/{T^ — T_) versus the 
number of particles in the FPU-/3 model for fixed and free boundary conditions. In the inset, 
one can see that the effective growth rate oje// defined in (106) is basically the same in both 
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Fig. 22. Thermal conductivity of the FPU-/? model versus lattice length N for r+ = 0.11, r_ = 0.09, 
and 6 = 1. The inset shows the effective growth rate a^ff versus N . Circles and diamonds correspond 
to free and fixed b.c, respectively. 
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Fig. 23. Thermal conductivity of the FPU-/3 model, r+ = 150, T_ = 15, fixed boundary condition. 
The data are taken from Ref. [76] 

cases, despite the clear diflFerences in the actual values of the flux itself. Additionally, a^ff does 
not deviate significantly from the theoretical prediction {a — 0.4). 

It is instructive to notice also that the now widely confirmed divergence of the thermal conduc- 
tivity with the chain length was already observed in previous simulations in spite of opposite 
claims made by the authors themselves. We refer to a paper by Kaburaki and Machida, where it 
was conjectured a slow convergence [76] of k{N) towards a finite value in the FPU-/3 model. By 
re-plotting their data in doubly logarithmic scales, a convincing power-law behavior is clearly 
seen instead, with even a quantitative agreement for the divergence exponent (see Fig. 23). 

Once the divergence is clearly established, the next question concerns the universality of the 
divergence rate. The discussion of this point involves considering a possible dependence on 

the temperature as well as on the leading nonlinearities [63]. Both questions arc addressed in 
Fig. 24, where n{N) is computed in the FPU-a model at a relatively low temperature. The 
convergence of the effective exponent towards 0.4 (see the inset in Fig. 24) suggests that the 
presence of a quadratic nonlinearity in the force field does not modify the overall scenario 
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Fig. 24. Thermal conductivity of the FPU-a model versus lattice length for g = 0.25, T = 0.1, and 
= 1. Triangles and circles refer to AT = 0.1 and AT = 0.02, respectively. The solid line corresponds 
to the linear divergence observed in a harmonic chain with the same temperatures. The inset shows 
the effective divergence rate CKg/ versus N for the data corresponding to full circles. 

observed in the FPU-/3 model. Additionally, notice that changes in the temperature gradient, 
without modifying the average T — (T^ + T_)/2, modify the effective conductivity only at 
relatively small sizes. In fact, we see in Fig. 24 that the two sets of measures corresponding to 
AT = 0.1 and 0.02 (triangles and circles, respectively) approach each other for larger than 
10^. In both cases k{N) increases linearly with N ior N < 10^ and no sizeable temperature 
gradient forms along the chain. Both facts hint at a weakness of anharmonic effects up to this 
time/length scales. This is confirmed by the comparison with the results for a pure harmonic 
chain (with the same setup and same parameters) that exhibit a clean linear growth of n with 
(see the solid line in Fig. 24) and a few-percent differences in the initial size range. The fact 
that K is smaller for larger AT can be thus attributed to a stronger boundary scattering that 
reduces the conductivity. 

The overall scenario is confirmed by the computation of k through the Green-Kubo formula. 
The very existence of an anomalous transport coefficient can be inferred from the slow decay 
of the corresponding autocorrelation function. In particular, if (J(i) J(0)) decays as t"^ with 
P < I, the integral in Eq. (118) diverges, thus signalling an infinite conductivity. Obviously, in 
finite chains one expects that an exponential decay eventually sets in, so that simulations should 
be performed for different chain lengths in order to be sure to pick up the truly asymptotic 
scaling behavior. From a numerical point of view, a simpler way to proceed consists in looking at 
the low- frequency divergence of the power spectrum S{u;) of the total heat flux: by the Wiener- 
Khinchin theorem, a power-law decay of the autocorrelations with exponent (3 translates into 
an o;^"^ behavior of S{u;) at small lo. 

The spectrum plotted in Fig. 25 is asymptotic in N in the selected frequency range. The low- 
frequency divergence of S{u;) imphes that the autocorrelation of J decays with a rate P ~ 0.63, 
i.e. that its time integral diverges. 

A quantitative comparison with the previous results can be performed by noticing that energy 
propagates with the constant sound velocity Vg- This can be understood by, e.g., looking at the 
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Fig. 25. Power spectrum S{uj) (in arbitrary units) of the global flux for an FPU-/? chain of length 
N = 1024 at a temperature T = 11.07. The curve results from an average over 1400 independent 
initial conditions. A blow-up of the low-frequency region is reported in the inset: the dashed line is a 
shifted fit with slope -0.37. 




Fig. 26. The spatio-temporal correlation function C{i,t) = (ii(i)io(O)) of the local flux for the FPU/3 
model. Micro-canonical simulations, energy density 8.8 

spatio-temporal correlation function C{i,t) = {ji{t)jo{0)) of the local heat flux [72] plotted in 
Fig. 26. Accordingly, one can turn the time divergence of k, as determined from the Green-Kubo 
formula into a divergence with N by restricting the integral in formula (118) to times smaller 
than the "transit time" Na/vg. This amounts to ignoring all the contributions from sites at a 
distance larger than N. With the above estimate of Cj, one obtains that k, oc N^'^. The latter 
exponent is the one reported in the last column of Table I. 

It is instructive to repeat the modal analysis for the contribution to the heat flux also in the 
nonlinear case. In spite of the fact that there are no longer eigenmodes, one can determine the 
contribution of each mode to the heat flux (the modes to be considered being defined according 
to the imposed boundary conditions) [87]. The relevant difference with the harmonic case is 
that, because of the lack of integrability, the A;-th mode does not only exchange energy with 
the left and right reservoirs, but also with all other modes (this is precisely the mechanism that 
is eventually responsible - in 3d - for a normal conductivity). Accordingly, one can write an 
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Fig. 27. Modal fluxes in an FPU chain of length 512 for fixed (solid line) and free (dashed line) 
boundary conditions. Vertical units are fixed in such a way that the total flux is normalized to unity 
in both cases. 

energy balance equation for the k-th mode as 



where is a self-explanatory notation for the fluxes towards the two heat baths, while J^' is 
the energy exchanged with the other modes as a consequence of the non-intcgrabic dynamics. 
Obvious global constraints imply that J2k ^fe ' = 0- I^irect simulations performed with both free 
and fixed boundary conditions suggest a much stronger property, namely that each "nonlinear" 
flux vanishes, JJ!^ = 0. The same simulations indicate that the effect of the boundary conditions 
on the modal fluxes is qualitatively similar to that in harmonic systems. In fact, from Fig. 27, 
we can see that for fixed boundary conditions, the contribution of low-A; modes is depleted 
and goes to for A; ^ 0, while a growth, if not a divergence, is observed for free boundary 
conditions. 

Prom the hydrodynamic description put forward in the previous chapter, one understands 
that the contribution to the anomalous behavior of heat conductivity arises from the behavior 
of long-wavelength modes in a similar way to the anomaly observed in disordered harmonic 
chains. It is therefore natural to ask why boundary conditions are so important in disordered 
harmonic chains that they may turn a diverging into a vanishing conductivity, while the same 
does not occur in nonlinear systems. The question becomes even more intriguing after having 
noticed that nonlinear systems are characterized by a similar dependence of the modal fluxes 
to that found in harmonic chains. A precise answer to this question would require combining in 
a single model internal relaxations (self-consistently described by mode-coupling theory) and 
dissipations due to the couphng with the external heat baths. Such a type of description is still 
lacking. 

Finally, we discuss the role of the boundary resistance in connection with the temperature 
dependence of conductivity. In fact, an interesting application of Eq. (56) has been proposed 
in Ref. [50] with reference to the FPU-/3 model. There, it has been empirically found that the 
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bulk conductivity scales with L and T as 



/ 1.2 L"T-i (T<0.1) , . 

''-\2L"rV4 (r>50) ■ ^ > 

According to kinetic theory, the conductivity can also be expressed as k — iVgCy. Since 
and Vs are almost constant and of order 1 in a wide temperature range, £ ~ k. Hence, at 
low temperatures the boundary jumps dominate the thermal profile up to the size that 
can be estimated according to Eq. (56). At low temperatures this effect is very strong since 
~ (2£/T)i^, while smaller boundary resistances are found at large temperatures, where 

~ (2£rV4)T^. 



6.3 The hard-point gas 



Although this review is basically devoted to analysing the behaviour of low- dimensional lattices, 
it is worth considering also fluid systems in so far as no qualitative differences are expected 
for the scaling behaviour of their transport properties. More specifically, in this section we 
discuss a set of point particles labelled by the index i — 1, ...N moving along a one- dimensional 
box extending from x = to x = L. The mass, position and velocity of the ith particle are 
denoted by rrii, Xi, and Ui, respectively. Interaction occurs only through elastic collisions. After 
a collision between the ith and the i + 1st particle, the respective velocities acquire the values 



mi-nii+i 2mi+i , 2mi rrii - rrii+i 

■ Ui^ ■ Ui+i , xij+i = ■ Ui ■ Ui+i, (148) 

rrii + rrii+i mi + rrii+i ^ mi + rrii+i rui + rrii+i 



as implied by momentum and energy conservation. Between collisions, the particles travel freely 
with constant velocity. Notice also that they mantain their initial ordering (no crossing is 
allowed) . 

The model is particularly suitable for numerical computation as it does not require integration of 
nonlinear differential equations. Indeed, the dynamics amounts simply to evaluating successive 
collision times and updating the velocities according to Eqs. (148). The only numerical errors 
are those due to round-off. The coupling with heat baths at the boundaries can be implemented 
in the usual way, e.g. by using Maxwell thermostats. Thus whenever a particle of mass m collides 
with a wall at temperature T, it is reflected back with a velocity chosen from the distribution 
P{u) = {m\u\/T)exp[-mu'^/{2T)]. 

In the limit where all the masses are equal, the system becomes integrable and one expects 
the same behavior observed in harmonic chains (see Sec. 4.1). In particular, the temperature 
proflle is flat (with T{x) = \/J\TZ), the heat current is independent of the system size and 
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no local equilibrium is attained. ^'^ However, as soon as the masses are different, the system 
is nonintegrable and (hopefully) ergodic, thus becoming a possible candidate for checking the 
validity of Fourier's law. Casati [88] considered the case of alternate masses {m2i = 1, m2i+i = 
(1 + 5) in suitable units ^''^ ). While from his numerical results it was not possible to draw any 
definite conclusion about the scaling behaviour of the conductivity, more recent simulations 
by Hatano [85] suggest a divergence rate consistent with what found for FPU (see Table I). 
Although we do not see any reason why coupled rotors and FPU should belong to different 
universality classes, the behaviour of conductivity in the hard-point gas appears to be still 
quite a controversial issue: simulations performed by Dhar [89] point to a slow divergence 
(k, ~ L°' with a < 0.2); equihbrium simulations discussed in Ref. [90] have led the authors 
even to conjecture a normal behavior; finally, further direct numerical studies confirm instead 
the existence of a divergence [92]. Accordingly, it seems reasonable to hypothesize that the 
uncertainty is due to strong finite-size effects that slow down the convergence to the expected 
asymptotic behaviour. 



Besides the dependence of k. on L, in Ref. [89] it has been studied the shape of the temperature 
profile. The simulations have been performed for different values of 6 and for up to 1281 
(adjusting the system size L so as to keep the average density of particles equal to 2). The 
number of particles has been chosen to be odd, so that the two particles in contact with the 
heat baths have always the same mass. Moreover, the averaging of the various observables 
has been performed over a time span corresponding to 10^ — 10^° collisions. In Fig. (28), it is 
plotted the steady-state temperature profile for different values of A^ and the same S = 0.22. 
The profile is smooth except at the boundaries, where two temperature drops are observed 
whose amplitude decreases with the system size. 



Upon increasing A^, the temperature profile approaches a limiting form. Quite amazingly, this 
shape is quite close to the one that would be predicted by kinetic theory. In fact, let us recall 
that kinetic theory applied to a one-dimensional gas predicts the Fourier law with a conductivity 
K, ~ VT. By then integrating the equation y/TdT/ dx — c.nt with suitable boundary conditions, 
one obtains 



Tk{x) = 



r+/'(i 



2/3 



(149) 



This corresponds to the solid curve in Fig. (28). The agreement with the numerical simulations 
is surprising, since kinetic theory predicts a finite conductivity. We are inclined to interpret this 
result as a confirmation of the existence of strong finite-size corrections. 



" In this respect, it is worth mentioning that the equal-mass case with dissipation has been studied 
by Du et al [91] as a toy model for a granular gas. They obtained a rather surprising stationary state 
which impUed a breakdown of usual hydrodynamics. 

From the invariancc of tlic dynamics described by Eq. (148) under mass rescaUng (rrij i^rrii), it 
follows indeed that the only independent parameters are the mass ratio (1 + 5) and the ratio of the 
boundary temperatures T^/T^. In fact the temperature profile does not change under rrii vrrii. 
Also from the boundary conditions, it is easily shown that T{yTj^, z/T_, x) = vT{T^, r_, x). 
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Fig. 28. Temperature profiles of the hard point gas for S = 0.22, r+ = 8,T_ = 2 and sizes 
= 21,41,81, 161,321,641 and 1281 (from top to bottom). The sohd Hne corresponds to Eq. (149). 
In the inset, it is plotted g{x) (see Eq. (150)) with the data for N = 161,321,641 and 1281 

The hard-point gas is interesting also for the possibility to investigate the nearly intcgrable 
regime when 6 <^ 1. The smaller is S, the larger has to be the system size in order to generate 
the same temperature profile. In Ref. [89] it has been conjectured that T{x, N, 5) depends on 
S and only through the scahng combination S'^N and it has been proposed the following 
scaling form 

T(x, N, 5) = n(x) + j^9{^)- (150) 

The above relation can be tested with reference to the data reported in Fig. (28). The good data 
collapse that can be appreciated in the inset of Fig. (28) supports the vahdity of the scahng 
relation (150) with 7 = 0.67. 

6.4 The coupled-rotor model 

The simplest example of a classical-spin Id model with nearest neighbor interactions lies in the 
class (2) with [93,94] 

l/(x) = 1 -cosx. (151) 

This model can be read also as a chain of N coupled pendula, where the pj's and the g^'s 
represent action-angle variables, respectively. It has been extensively studied [95,96,97] as an 
example of a chaotic dynamical system that becomes integrable both in the small and high 
energy limits, when it reduces to a harmonic chain and free rotors, respectively. In the two inte- 
grable limits, the relaxation to equilibrium slows down very rapidly for most of the observables 
of thermodynamic interest (e.g., the specific heat) [96,97]. As a consequence, the equivalence 
between ensemble and time averages is established over accessible time scales only inside a lim- 
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Fig. 29. Conductivity k versus chain length N as obtained from non-equihbrium molecular dynamics. 
Circles correspond to the rotator model with temperatures T+ = 0.55 and T_ = 0.35. The dashed line 
represents the best fit with the function a + b/N. The shaded region represents the uncertainty about 
the conductivity on the basis of the Green-Kubo formula. 

ited interval of the energy density e. Here, we focus our attention mainly on heat conduction 
in the strongly chaotic regime. 

In Refs. [93,94], it has been shown that, contrary to the expectations, this model exhibits a 
finite conductivity in spite of the existence of an acoustic branch in its spectrum in the harmonic 
limit. In Ref. [93], simulations have been performed for T+ = 0.55, TL = 0.35, and chain lengths 
ranging from A?" = 32 to 1024 with fixed boundary conditions and Nose-Hoover thermostats. 
The equations of motion have been integrated with a 4-th order Runge-Kutta algorithm and 
a time step At = 0.01. The results, reported in Fig. 29 clearly reveal a convergence to a value 
of K approximately equal to 7 (see the circles). The dotted line in good agreement with the 
numerical data is the best fit with the function a + b/N^'^ . However, more important than 
assessing the convergence properties of k,{N) is to notice its finiteness for N ^ oo. 

In fact, this is the first system where normal heat conduction has been convincingly ascertained 
in the absence of an external field. Precisely because of this atypical behavior, it is important 
to confirm this result with a computation of thermal conductivity through the Green-Kubo 
formula. In Ref. [93] , micro-canonical simulations have been performed in a chain with periodic 
boundary conditions. In the absence of thermal baths, the equations of motion are symplectic; 
accordingly, the Authors have made use of a 6-th order McLachlan-Atela integration scheme 
[33], fixing the energy density equal to e = 0.5, a value that corresponds to T pa 0.46, close 
enough to the average temperature in the nonequilibrium simulations. 

The correlation function has been computed by exploiting the Wiener-Khinchin theorem, i.e. 
by anti-transforming the Fourier power spectrum of the total flux. The curve plotted in Fig. 30 
indicate both a clean exponential decay and an independence of the behavior of for N > 256 
(at least in the reported time range). This allows for an accurate determination of the integral 
of Cj{t). The gray region in Fig. 29 corresponds to the resulting value of k, taking into account 

Notice that this function is asymptotically equivalent to expression (56) , derived under the assump- 
tion that finite-size corrections are due only to the boundary resistance. 
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Fig. 30. The autocorrelation function of the total heat flux in a chain of coupled rotors with periodic 

boundary conditions and energy density e = 0.5. Dashed, dot-dashed and dotted lines correspond to 
N = 256, 512, and 1024, respectively. The solid line, corresponding to Cj = exp(— i/30) has been 
drawn for reference. 

statistical fluctuations. The quantitative agreement between the two estimates of the heat 
conductivity is important in that it confirms also the finiteness of k in a context where this was 
not a priori obvious. 

In order to emphasize the difference between the dynamics of the present model and that of the 
previous systems, it is instructive to look at the power spectrum of the low-A; Fourier modes. 
In Fig. 31 it is possible to compare the spectra of some low-k modes in coupled rotors with 
those in a diatomic FPU-/3 chain. In the latter case, sharp peaks are clearly visible (notice 
also that the peaks become increasingly narrow upon decreasing k): this is a signal of an 
effective propagation of correlations [72]. Conversely, in the rotors, the low-frequency part of 
the spectrum is described very well by a Lorcntzian with half- width 7 = Dk^ (D ^ 4.3). This 
represents an independent proof that energy diffuses, as one expects whenever the Fourier's law 
is estabhshed. 

In the attempt to explain the striking difference in the transport behavior exhibited by this 
model with respect to that of the previous models in the same class, one cannot avoid noticing 
that the pair potential V{qi^i — qi) possesses infinitely many equivalent valleys. As long as 
{qi+i — qi) remains confined to the same valley, there is no reason to expect any qualitative 
difference with, e.g., the FPU-/9 model. Phase slips (jumps of the energy barrier), however, 
may very well act as localized random kicks, that contribute to scattering of the low-frequency 
modes, thus leading to a finite conductivity. In order to test the validity of this conjecture, one 
can study the temperature dependence of k, for low temperatures when jumps across barriers 
become increasingly rare. The data plotted in Fig. 32 indicate that the thermal conductivity 
behaves as k ~ exp{r)/T) with r) ~ 1.2. The same scaling behavior is exhibited by the average 
escape time r (see triangles in Fig. 32) though with a different 77 ^ 2. The latter behavior can be 
explained by assuming that the phase slips are the results of activation processes. Accordingly, 
the probability of their occurrence is proportional to exp(— AV"/T), where AV is the barrier 
height to be overcome. The behavior of r is thus understood, once we notice that AV = 2. In 
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Fig. 31. Power spectra of the 1st, 2nd, 4th, 8th, 16th, 32nd, and 64th Fourier mode in arbitrary units 
for a chain oi N = 1024 particles. Panel (a) refers to a chain of rotors with energy density e = 0.5 
(the wavenumber increases from left to right); panel (b) refers to a diatomic FPU-/? chain with masses 
1, 2 and energy density e = 8.8 (the wavenumber increases from top to bottom in the low-frequency 
region. In both cases the curves result from an average over 1000 independent simulations. 




Fig. 32. Thermal conductivity k versus the inverse temperature 1/T in the rotor model (open circles). 
Triangles correspond to the average time separation between consecutive phase slips in the same 
system. 



the absence of phase slips, the dependence of the conductivity on the length should be the same 
as in FPU-systems, i.e. k, N'^^^. In the presence of phase slips, it is natural to expect that the 
conductivity is limited by the average distance between consecutive phase slips. Under the 
further assumption of a uniform distribution of the slips, their spatial and temporal separation 
has to be of the same order, thus implying that k(T) exhibits the same divergence as r for 
T — i> 0, though with a different rate k, ~ exp[2AI//(5T)]. Therefore, at least on a qualitative 
level, one can indirectly confirm that phase slips are responsible for the normal heat transport. 
On a quantitative level, however, there is a discrepancy between the observed and the expected 
value of the exponent t] (1.2 vs. 0.8). Among the possible explanations for the difference, we 
mention the presence of space-time correlations in the pattern of phase-slips and the existence 
of ever increasing deviations from the asymptotic law k, ~ N"^/^ for T — > 0, due to the vanishing 
of nonlinear terms. 
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In order to further test the conjecture that jumps between adjacent vaUeys of the potential are 
truly responsible for a normal heat transport, some other models have been studied as well. 

Let us start by considering an asymmetric version of the rotor model, namely 

V{x) ^A-cosx + 0.4 sin 2x (152) 



where A is fixed in such a way that the minimum of the potential energy is zero. Simulations 

performed at a temperature corresponding to one quarter of the barrier-height again indicate 
that the conductivity is finite, confirming the empirical idea that jumps are responsible for 
breaking the coherence of the energy flux [93] . 

It is instructive to look more closely at the behavior of this model. In view of the asymmetric 

potential, one might expect that the average force (j) = {J2fi)/N is non-zero (like, e.g., in the 
FPU-a case). Nonetheless, micro-canonical simulations show that although the distribution of 
forces is definitely asymmetric, their average value is numerically zero. This can be understood 
by noticing that in view of the boundedness of the potential, the system cannot withstand any 
compression. 

Accordingly, we are led to introduce yet another model where the existence of more than 
one valley is accompanied by an unbounded potential. The simplest way to achieve this is by 
considering the double- well potential 

V{x) = -xV2 + xV4 (153) 



(it the same as in FPU-/3 with the opposite sign for the harmonic term). The results of direct 

simulations [98] performed with Nose-Hoover thermostats are reported in Fig. 33 for three 
different values of the temperature below the barrier height. One can see that the growth of the 
conductivity, after an initial slowing down, increases towards, presumably, the same asymptotic 
behavior observed in the FPU model. This is at variance with the preliminary simulations 
reported in Ref. [93], where it was conjectured a normal heat transport. These results thus 
indicate that jumps alone from one to another valley are not sufficient to destroy the coherence 
of low- A; modes dynamics. It is necessary that subsequent jumps be independent of one another. 
In the double-well potential (153), gj+i — Qi cannot exhibit two consecutive jumps to the right. 



7 Anharmonic chains with external substrate potentials 



In this section we consider the class of models described by the Hamiltonian (5), with a non- 
vanishing substrate potential U. This means that translational invariance breaks down and 
momentum is no longer a constant of the motion. Accordingly, the dispersion relation is such 
that uj{q)j^O for g = 0. 
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Fig. 33. Thermal conductivity k versus chain length for the potential (153) and three different av- 
erage temperatures: 1/16 (circles), 3/32 (squares), and 1/8 (diamonds). In all cases the temperature 
difference is chosen equal to 1 tenth of the average temperature. 

7.1 Ding-a-ling and related models 

The so-called ding-a-ling model was first introduced by Dawson [99] as a toy model for a 
Id plasma. It can refer to different contexts: (i) a set of identical charge-sheets embedded in a 
fixed neutralizing background; (ii) a system of harmonic oscillators with the same frequency and 
equilibrium positions sitting on a periodic lattice and undergoing clastic collisions that exchange 
their velocities. Notice that in the low-energy limit, it reduces to the Id Einstein crystal, i.e. 
set of independent harmonic oscillators all having the same frequency (no dispersion) . 

Independently of [99], Casati ct. al [100] introduced a modified version, where the harmonic os- 
cillators (say the even-numbered particles) alternate with free particles of the same (unit) mass. 
The latter are only constrained to lie between the two adjacent oscillators. The Hamiltonian 
can be symbolically written as 



where ui — u for even / and zero otherwise. A common feature of this class of models is that 
within collisions the motion of the particles can be determined analytically so that the basic 
requirement is the computation of the occurrence times of the collision events. Therefore, the 
dynamics naturally reduces to a discrete mapping. 

For the isolated system (e.g. a chain with periodic boundary conditions) the dynamics depends 
only on the dimensionless parameter e = e/iujaY where e is the energy per particle and a 
the lattice spacing. The Authors of Ref. [100] studied the dynamical behavior of the model by 
fixing e — 1 and changing uj. They concluded that, for uj and N large enough, the dynamics is 
strongly chaotic and soliton-like pulses are sufficiently attenuated [88]. This renders the model 
a good candidate to check the validity of Fourier's law. 

Finite thermal conductivity - The validity of the Fourier's law was first established by perform- 
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ing a scries of non-equilibrium simulations, where the freely moving end-particles were put in 
contact with two Maxwellian reservoirs. The average flux J was then computed by summing 
the amounts of energy 5E exchanged with one of the reservoirs in all collisions during the 
simulation time. The average temperature gradient was estimated with a linear fit (to get rid 
of boundary effects) . By evaluating the thermal conductivity as a function of the lattice length 
up to A/" = 18 for T+ = 2.5, T_ = 1.5 and a; = 1, it was concluded that k{N) attains a constant 
limiting value already ioY N > 10 . 

After having estabhshed the existence of a finite value of the transport coefficient, the Authors 
have compared the value of k with the result of linear response theory. To this aim. because of 
the discontinuities due to the collision processes, it was preferred to express the Green-Kubo 
formula in terms of the integral quantity 

t+T 

AQ{t,T)^ j J{to)dto 
t 

From Eq. (118), recalling that the ensemble average is equivalent to a time average (and under- 
standing the limit N —>■ oo), it is straightforward to show that (for more details see Ref. [103]) 

^OKiT) = lim^^{{AQit,r)r\ , (155) 

where the subscript t indicates that the average is performed over the time variable t. 

Additionally, it was decided to compute the total heat fiux not by summing up the pf local 
contributions as, for instance, done in Ref. [90] for the hard point gas, but, more directly, 
determining AQ{t, r) as the amount of energy exchanged in all collisions occurred in the interval 

[t,t + T]. 

Casati et al. provided a convincing numerical evidence that the limit (155) exists for a closed 
chain of 48 particles and cu = 10. In this case, the energy transport is diffusive and they showed 
that the Kci^-value obtained in this way is in good agreement with the one obtained from direct 
simulations. 

The results of Casati et al. have been lately reconsidered by Mimnagh and Ballentine [101] 
who performed a detailed series of simulations with longer chains and in a wider parameter 
range. Curiously enough, they found that the value of k reported in Ref. [100] is not the true 
asymptotic value (achieved only for N > 200) but a minimum of k{N). This can be seen in 
Fig. 34, where one can also appreciate how the unfortunate choice of working with N < 20 may 
give the false impression of a convergence towards a smaller value. This fact does not affect 
the correctness of the conclusion reported in Ref. [100] and the importance of the results, but 

Recall, indeed, that the substrate potential does not contribute to the fiiix. 
^'^ As for the similar comparisons discussed in the previous chapter with reference to the coupled 
rotors, micro-canonical simulations have to be performed for an energy density that corresponds to 
the average temperature in the non-equilibrium simulations. 



68 



50 



30 



20 



10, 





50 

r ^ 

1 / 


















/ 30 
" 20 

10 




















10 


15 N 20 



200 



400 



600 



Fig. 34. Thermal conductivity of the ding-a-hng model. Size dependence of «; for a; = 1 and e = 1.5 
(from [101]). In the inset, it is presented an expanded view in the range of sizes considered in Ref. [100] 
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Fig. 35. The coefficient /x^ as a function of e for the ding-a-ling model (from [101]) 

signals again how cautious one should be in drawing conclusions from the study of relatively 
short systems! Motivated by this observation, Mimnagh and Ballcntinc carefully studied finite- 
size corrections in a wide range of £- values, by plotting the resistivity p(A^) = 1/ k{N) versus 
N . In all cases, the data are well fitted by 



p{N) = Poo 1 + 



(156) 



for N large enough. Accordingly, the minimal chain length required to obtain an estimate of 
Poo with a fixed relative accuracy is proportional to p^, a quantity which is found to increase 
dramatically for e > Ec = 0.04 (see Fig. 35). This is readily understood as e — > cxo is an 
integrable limit that corresponds to a gas of bouncing free particles. The asymptotic resistivity 
Poo is a monotonously decreasing function of e, displaying a crossover from a slower to a faster 
decay at Ec (see Fig. 36). Since a similar crossover is found when looking at both the collision 
rate and the maximum Lyapunov exponent, it is natural to interpret the phenomenon as a sort 
of transition from strong to weak chaos, in close analogy with what found, e.g., in the FPU 
model [102]. 

Finally, as for the temperature profile, the nonlinear shape reported in Fig. 37 reveals a sizeable 
temperature dependence of the conductivity. 

Modified models - Prosen and Robnik [103] considered the original Dawson system where, 
as already mentioned, the free particles are removed (i.e. they set uoi = 1 for all particles in 
(154)). They named it ding-dong model to distinguish it from the one presented above. Their 
careful numerical study confirmed the validity of Fourier's law in a wide temperature range. 
Besides direct non-equilibrium simulations with Maxwellian thermostats and the Green-Kubo 
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Fig. 36. The asymptotic resistivity Poo versus e (from [101]) 
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Fig. 37. Temperature profile for the ding-a-hng model, uj = 2 (from [101]) 



formula, they also implemented an efficient transient method that allowed them to explore the 
high temperature regime (T > 3), where, because of the nearly integrable dynamics, a slow 
convergence of the averages with time and/or size is observed. In the opposite, low-temperature, 
limit (T < 0.1), they are able to prove that the conductivity vanishes as exp(— 1/47") finding a 
reasonable agreement with numerical results. 

Posch and Hoover [104] investigated a modified ding-a-ling model where the harmonic potential 
is replaced by a gravitational one 



where gi is a constant acceleration for even / and zero otherwise. Their simulations further 
confirm that thermal conductivity is finite for this class of models. Moreover, they computed 
the spectrum of Lyapunov exponents in the non-equilibrium steady state, showing that the 

microscopic dynamics takes place on a strange attractor. Finally they observed that the heat 
flux is proportional to the difference between the phase-space dimension and that of the strange 
attractor. 

Finally, it is worth mentioning the exactly solvable model studied by Kipnis et al. [105] that 
could be regarded as the "stochastic version" of the ding-ling model. As in the latter, it consists 
of a linear array of harmonic oscillators but the interaction occurs via a random redistribution 
of the energy between nearest neighbors rather than through deterministic collisions. More 
precisely, upon denoting with = + qf the energy of the l-th oscillator (in suitable units). 



7i = — \- mgi\qi\ + "hard point core' 




(157) 
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the dynamics is given by the updating rules 

ei = ^[+1 = , (158) 

where P is a random variable uniformly distributed in the interval [0, 1]. The total energy is 
thus kept constant except for the two oscillators at the chain extrema, that are in contact 
with reservoirs at different temperatures T± according to a Glauber dynamics. The Authors 
were able to rigorously show that a unique stationary non-equilibrium measure exists and to 
compute both the temperature profile and the heat flux in the steady state: 

T{x) = T_ (^^) + T+ (^-^) , -l<x<l (159) 



j = -^(r^-r_) . (160) 

The last results imply that Fourier's law holds and that the thermal conductivity is equal to 
A;b/2 in the chosen units. 



7.2 Klein- Gordon chains 



An important subclass of models (5) is the one in which the inter-particle potential is harmonic 

^ = E [|^' + U{qi) + Iciqi^i - qir] ; (161) 

it is often referred to as the Klein-Gordon lattice. The latter has recently received a great 
attention as a prototypical system where strong discreteness effects may come into play in the 
limit of small C. 

The first and most complete study of the transport problem in this class of models has been 
carried on by Gillan and HoUoway [42] for the Frenkel-Kontorova potential 

U{x) = -Uo cos (^^] . (162) 



a 



The model can be interpreted as a chain of either coupled particles in an external periodic field 
or torsion pendula subject to gravity. In the latter case a = 27r and qi represents the angle with 
respect to the vertical direction: it can be read as the discretized (and non-integrable) version 
of the well-known sine- Gordon field equation. 

Besides energy, the dynamics admits a further conserved quantity, the winding number V, which 
is an integer defined by the boundary condition g^+jv = qi + dV. In the particle interpretation V 
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Fig. 38. Thermal conductivity versus chain length in 0'^ chains with Nose-Hoover thermostats (0 = 10). 
Panel (a) refers to the single well case (a = 6 = 1 in Eq. 163): the results have been obtained for C = 1., 
r+ = 8, and T_ = 6. The shaded region represents the value obtained from the Green-Kubo formula 
with its statistical uncertainty. Panel (b) refers to the double-well case (a = —1, b = 1) for an average 
temperature T = 0.37 and a temperature difference 0.002. The dashed line is just a guide for the eyes. 

represents the number of potential wells, while for the pendula it can be viewed as the degree 
of built-in twist in the system. According to the general theory of irreversible processes [1], 
transport involves thus the flow of both particles and energy caused by gradients of the number 
density and temperature. Accordingly, a 2 x 2 matrix of transport coefficients is required, but 
because of the Onsager relations, there exist only 3 independent transport coefficients, that 
are chosen to be the thermal conductivity, the diffusion coefficient and the heat of transport. 
Gillan and HoUoway computed the thermal conductivity numerically in the general case of 
non-vanishing winding number with three different methods: (i) attaching two heat baths; (ii) 
through the Green-Kubo formula; (iii) by adding an external field (see Chap. 3). All the methods 
give consistent results and clearly indicate that the thermal conductivity is finite. 

Their results were later confirmed by a numerical study by Hu et al. [106] who investigated the 
dependence of the transport coefficient on the lattice length (for P = 0). In Ref. [106] it was 
also shown that the same holds for a more general version of the Frenkel-Kontorova model with 
an anharmonic inter-site potential. 

In order to illustrate the type of behavior that is observed in this class of systems we show in 
Fig. 38 some data for the 0*^ chain 



In panel (a), we present a case of fast convergence to a small k value for a single- well po- 
tential; panel (b) refers instead to a low-temperature regime characterized by large thermal 
conductivity. 

In physical terms, one has to apply equal and opposite forces (in the particle version) or torques 
(pendulum version) to the two end particles, in order to maintain the required value of V which can 
be also interpreted as the net number of kink excitations (i.e. the number of kinks minus the number 
of anti-kinks). 



U{x) 




(163) 
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Evidence of a finite conductivity for the case a — has been reported in Refs. [107,49]. The 
conductivity is found to decrease with temperature according to the law [49] 



(164) 



that is reminiscent of what often experimentally observed in insulating crystals. Upon increasing 
AT, there exists a "transition" to a nonlinear regime characterized by a non uniform local 
temperature gradient. By including the empirical law (164) into Fourier's law, it is possible to 
check the consistency with the measured T{x): a good agreement with the simulations is found, 
provided boundary jumps are taken into account. 

Finally, Tsironis et al. [108] obtained further numerical evidence of the existence of a finite 
thermal conductivity for systems like (161). Beside reconsidering the Frenkel-Kontorova po- 
tential (162), two further examples were analyzed, the sinh-Gordon and bounded single- well 
potentials. 



as representatives of the classes of hard and soft anharmonicity, respectively. 



8 Integrability and ballistic transport 

When the equilibrium dynamics of a lattice can be decomposed into that of independent 
"modes", the system is expected to behave as an ideal conductor. The simplest such example 
is obviously the harmonic crystal, that has infinite conductivity and cannot, therefore, sup- 
port any temperature gradient. However, this applies also to the broader context of integrable 
nonlinear systems. They are mostly one-dimensional models characterized by the presence of 
"mathematical solitons" , whose stability is determined by the interplay of dispersion and non- 
linearity. This interplay is expressed by the existence of a macroscopic number of conservation 
laws constraining the dynamical evolution. Thereby, the existence of stable nonlinear excita- 
tions in integrable systems is expected to lead to ballistic rather than to diffusive transport. As 
pointed out by Toda [109], solitons travel freely, no temperature gradient can be maintained 
and the conductivity is thus infinite. From the point of view of the Green-Kubo formula, this 
ideal conducting behavior is refiected by the existence of a nonzero fiux autocorrelation at ar- 
bitrarily large times. According to the discussion reported in Section 6.2, this, in turn, implies 
that the finite-size conductivity diverges linearly with the size. 

Although integrable models are, in principle, considered to be exactly solvable, the actual 
computation of dynamic correlations is technically involved. A more straightforward approach 
is nevertheless available to evaluate the asymptotic value of the current autocorrelation. This 
is accomplished by means of an inequality due to Mazur [110] that, for a generic observable A, 





(165) 



73 



reads as 

}^J-}{Ammdt>^^-^ , (166) 

where (...) denotes the (equihbrium) thermodynamic average, the sum is performed over a set 
of conserved and mutually orthogonal quantities Qn {{QnQm) — {Qn)^n,m)- Furthermore, it is 
assumed that (A) — 0. 

Zotos [111] applied the above result to the equal-masses Toda chain with periodic boundary 
conditions, defined, in reduced units, by the Hamiltonian 

(167) 



1=1 



pt 



+ exp(-r;) 



where r/ = xi+i — xi is the relative position of neighboring particles. As is known [112], the 
model is completely integrable as admits N independent constants of th motion, the first among 
which are 



N 



Qi^Y^Pi (168) 
1=1 

N 2 

Q2 = T.^i + e-^^ (169) 

N 3 

Q^ = Y.^-k+{Pi+Pi+i)e-'' (170) 

1=1 ^ 

^ 1 
^4 = E T + ^P'l + + ^^Hi)^"'' + 0^"'" + e-'''e-'-'+^ (171) 

N 5 

Q^ = Y.-k + {p'+ PiPi+1 + PiPli + (172) 

1=1 ^ 

+ iPi + Pi+i)e-^'' + {Pi + 2pz+i + pi+2)e-'-'e-'^+K.. (173) 

The "trivial" conserved quantities Qi (the total momentum) and Q2 (the total energy) are 
of course present in all translationally invariant systems of the form (2), irrespective of their 
integrability. 

Let us consider the fixed temperature-pressure thermodynamic ensemble, 

L 

{A{t)A{0)) ^Z-^ l[ dpidriA{t)A{0)e-^^'^+^''^ (174) 
1=1 



where Z — J Ylj^i dpidrie P(.'H+pl)^ ^ _ J2iLiri is the "volume" of the chain and P is the 
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pressure. In this ensemble, equal-time correlation functions can be calculated analytically. For 
instance, the average inter-particle distance is 

(r) = \n(3-^f((3P) (175) 

where ^(-z) is the digamma function. 

The total heat flux is given by (see Eqs. (23) and (30) ) 

I 

2 

where hi — ^ + |(e~''' + e"^'-^). In order to apply the Green-Kubo formula in the chosen 
ensemble, one has to consider a "shifted" flux (see the discussion at the end of Sec. 5.2) 

This is equivalent to removing the contribution of Qi in the right hand side of the Mazur 
inequahty (166) for ^4 = J. 

Lower bounds on the long time value of (J(t)J(O)) can thus be calculated by the inequality 
(166), using the first m conservation laws [111]. As has a structure very similar to the energy 
current, a large contribution from this term has to be expected. Moreover, Qn with even n are 
uncoupled with J, so that it suffices to consider odd values of n. 

In order to utilize the inequality (166), it is not necessary to orthogonalize the conserved 
quantities (168-173). One can, indeed, replace the sum of the first m terms in the r.h.s. of 
Eq. (166) with 

= {J\Q){Q\Q)-'{Q\J) (178) 

where {Q\Q) is the mxm overlap matrix of Qn'n and {Q\J) is the overlap vector of J with the 
Q'^s. The ratio C"*/ ( J^), representing a lower bound to the conductivity is found to increase 
monotonously with the temperature. At low T, the growth is linear with a slope comparable 
to the density of solitons Ns/N = ln2/7r^T. This trend is interpreted as an evidence for the 
increasing contribution of thermally excited nonlinear modes to ballistic transport. 



9 Two-dimensional lattices 

It is well known that many properties of statistical systems depend on the dimension d of the 
space, where they are embedded. In Chapter 5, we have seen that transport properties are not 



, , (Pl+i +Pl) -n 



(176) 
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expected to violate this rule. In this Chapter we discuss some results of molecular-dynamics 
simulations in two dimensions. In fact, as soon as the dimension of the physical-space is set to 
a value larger than 1, the direct investigation of sufficiently large systems becomes problematic. 

First, we briefly discuss the numerical studies appeared so far in the literature, presenting 
them in a historical perspective. Then, we more extensively discuss some recent numerical 
experiments that have allowed verifying the predictions of mode-coupling theory also in two 
dimensions. 



9.1 Early results 

To our knowledge, the first attempts of investigating the heat conduction problem in 2d lattices 
with (at that time) heavy numerical simulations are two papers by Payton, Rich and Visscher 
[32,113], that appeared more than three decades ago. These Authors investigated the combined 
effect of nonlinearity and disorder on heat conduction in 2d harmonic and anharmonic lattices, 
with Lennard- Jones pair potentials. Their studies aimed also at analyzing the dependence of 
heat conductivity k on the temperature and on the concentration of impurities, as a measure 
of disorder. They found evidence of an increase of k in disordered nonlinear systems compared 
to the harmonic case (see also Chap. 6). 

On the other hand, the dependence of k on the system size was ignored, probably because the 
Authors did not consider this a problem of major concern. In fact, according to the classical 
view of Peierls [2], phonon-phonon scattering processes were assumed to be sufficiently efficient 
to ensure normal transport properties in the presence of strong nonlinearity and disorder. 

Later, Mountain and MacDonald [114] performed a more careful study on the dependence of k 
on the temperature T. They considered a 2d triangular lattice of unit-mass atoms, interacting 
via a Lennard- Jones 6/12 potential. At variance with the previous investigations, no disorder 
was included, and their numerical results were consistent with the expected classical law k ~ 
T~^. Again, the dependence of k on the system size was not investigated. 

The first contribution in this direction is the paper by Jackson and Mistriotis [79]. These 
Authors compared measurements of k in the Id and 2d FPU lattices: they concluded that in 
both cases there was no evidence that the transport coefficient is finite in the thermodynamic 
limit. It is worth mentioning an interesting remark by these Authors: "the dependence of k on 
the system size cannot be adequately described in the high temperature, i.e., classical, limit by 
Peierls' model of the diluted phonon gas, because the perturbative Umklapp processes cannot 
account for the genuine nonlinear effects that characterize such a dependence" . 

Conversely, more recent molecular-dynamics simulations of the 2d Toda- lattice [115] have been 
interpreted in favor of the finiteness of k in the thermodynamic limit. Further confirmation in 
this direction can be found in appendix B of the interesting paper by Michalski [116], where 
heat conductivity in models of amorphous solids was thoroughly investigated. 
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The interplay of disorder and anharmonicity, that inspired the first contributions by Payton, 
Rich and Visscher has been reconsidered by Poetzsch and Bottger [ff7,f f8] who investigated 
percolating and compositionally disordered 2d systems. In particular, they tried to pinpoint 
the role of third- and fourth-order anharmonicity, concluding that, at equal temperature, the 
latter yields a larger value of k than the former one. Moreover, the Authors reported also about 
the dependence of k on the system size. They assumed the same point of view of Michalski 
[116], but a careful inspection of Fig. 5 in [117] shows that their data are also compatible with 
a systematic increase of k, with the system size. 

In a more recent and accurate investigation, Dellago and Posch [119] studied, by molecular- 
dynamics techniques, heat conduction as well as Lyapunov instability in a generalized version 
of the XF-model. Besides various interesting results, the manuscript contains a very clear 
indication of the finiteness of k in the thermodynamic limit. In the light of what discussed 
in Section 6.4, this result is not surprising, since k is finite already in Id for models of this 
type. A further interesting remark contained in Ref. [119] concerns the behavior of k below the 
transition temperature of the XF-model when the diffusive behavior of energy transport is lost 
and anomalous behavior seems to set in. 

9.2 Divergence of heat conductivity 

Heat conduction in 2d models of oscillators coupled through anharmonic, momentum-conserving 
interactions is expected to exhibit different properties from those of Id systems. In fact, MCT 
predicts a logarithmic divergence of k with the system size N at variance with the power-law 
predicted for the Id case (see Section 5.3). 

Following [69] , we discuss the results of molecular-dynamics simulations of the FPU-/3 potential 
(see Eq. (4), g^, = 0) and the LJ-(6/12) potential (see Eq. (3)). For the sake of simplicity we 
introduce the shorthand notations Vi{z) and V2{z) to denote, respectively, the two models. The 
goal of this twofold choice is to verify that the prediction of MCT is truly independent of the 
potential, provided it belongs to the class of anharmonic momentum-conserving interactions. 

While Vi does not contain natural scales for both distances and energies, the natural length scale 
of V2{z) is the equilibrium distance a, while its energy scale is the well depth e. Therefore, after 
having arbitrarily fixed g2 — i and g^ — 0.1 in Vi{z), a and e have been determined (a = 25, e8.6) 
by imposing that the coefficients of the second and fourth order terms of the Taylor expansion of 
V2 around its minimum coincide with g2 and g^, respectively. Notice, however, that, at variance 
with Vi{z), the Taylor expansion of V2{z) incorporates also a non- vanishing cubic term. Despite 
both models are characterized by the same shortest harmonic time scale, Tmin = 7r/\/2, the L-J 
model has been integrated with a slightly shorter time step (At = 5 • 10~^, compared to 10~^ 
for the FPU model) in order to ensure a sufficient accuracy when the dynamics experiences the 
strong nonlinearities of V2{z) (nearby the divergence of the Lennard- Jones potential). 

In both cases, simulations have been performed with reference to a square lattice containing 
Nx X A^j^ atoms of equal masses m. The equilibrium position of each atom has been chosen 
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to coincide with a lattice site, labelled by a pair of integer indices The origin of the 

Cartesian reference frame is fixed in such a way that 1 < i < A^^, and 1 < j < Ny. Accordingly, 
the components of the 2d- vector of the equilibrium position, r^j, are given by the interger pairs 
Moreover, in analogy to the Id case, each atom has been assumed to interact with its 
nearest-neighbors (herein identified with the von Neumann neighborhood). 



The model is thus represented by the Hamiltonian 

> |2 



n 



EE 

i=l 3=1 



I 



2m 



+ ^(iQi+ij - %-|) + - Qijl) 



(179) 



where qij(t) = rij(t) — r^j, i"ii(^) is the instantaneous position vector of the (i,j)-atom, and 
Py(i) is the corresponding momentum vector. 

For what concerns the definition of local temperature, the following definitions are all equivalent 
(see also Sec. 2.2), 



ly)\2 



,(a:)\2 



m 



m 



2m 



;i80) 



where p\^^ and p^J^ are the x and y components of the momentum vector pj^ , respectively. 

Moreover, tedious but straightforward calculations, akin to those presented in Section 2.3 , 
allow to express the x and y components of the local heat flux as 



■ (x) 
Jij 

Ay) ^ 

Jij 



Am 
a 

4m 



fxx ( Jx) (x) \ r-yx ( (y) iy) \ 

Jij [Pij + Pz+lj) + Jij [Pij + Pi+lj) 

fxy ( (x) , ix) \ ryy / (y) (y) \ 

Jij [Pij + Pij+1 ) + Jij [Pij + Pij+1 ) 



where the components of the local forces are defined as 

fxx ^ dy {\^i+lj - ^ij\) ryx _ _ 
Jij „ Ct> J 'I 



dV {\q,i+ij - Cii 



fxy ^ 9V {\q,ij+i - q,ij\) yy 

•lij o {x) Ji' 



dV {\q,ij+i - q,ij\) 

dqfj' dJif 
Finally, in analogy with Eq. (30), the total heat flux vector is deflned as 



J - EJ 



(181) 



(182) 



^,3 



The non-equilibrium simulations have been performed by coupling all atoms on the left (right) 
edge of the 2d lattice with the same thermal bath operating at temperature (T_). In the 
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Fig. 39. Heat conductivity k versus the system size for the 2d FPU (3 (a) and Lennard-Jones (b) 
models. In panel (a) r+ = 20 and T_ = 10; in panel (b) T^. = 1 and r_ = 0.5. In both cases, statistical 
errors have the size of the symbols. 

numerical studies reported hereafter, thermal baths have been simulated by applying the Nose- 
Hoover method. Nonetheless, the Authors of Ref. [69] verified that the same results are obtained 
upon using stochastic thermal baths, as well. Periodic and fixed boundary conditions have been 
adopted in the direction perpendicular (y) and parallel {x) to the thermal gradient, respectively. 
Finally, let us notice that, for the investigation of the thermodynamic limit, the simulations for 
different lattice sizes should be performed by keeping the ratio R = Ny/N^ constant. From the 
numerical point of view, it is convenient to choose small R values, since for a given longitudinal 
length aNx, the simulations are less time consuming. However, too small ratios would require 
considering larger system sizes to clearly observe 2d features. In [69] it was checked that R — 1/2 
is a good compromise for both Vi and V2 choices. 

With the above physical setup, heat equation (1) implies that a constant thermal gradient 
should establish through the lattice in the x-direction with (J^^^) > and (.1^^^) = 0. The time 
span needed for a good convergence of the time-average (■) increases with N^: for instance, 
C(IO^) units proved sufficient for — 16, while O(IO^) units are needed when — 128 (for 
not too small energy densities). 

The detailed analysis of temperature profiles performed in Ref. [69] has revealed deviations from 
the linear shape predicted by Fourier law (this is particularly true in the case of the Lennard- 
Jones potential V2), but one cannot exclude that this is to be attributed to the relatively large 
temperature differences adopted in order to have nonnegligible heat fluxes. Anyway, despite 
such deviations, simulations provide convincing evidence that the temperature gradient scales 
like A^""*^. Accordingly, the dependence of k, with the system size aN^ can be determined by 
plotting K oc {j^^^)Nx versus A^^;. The data reported in Fig. 39 support the MCT prediction of 
a logarithmic growth both for the FPU and Lennard-Jones potentials for two rather different 
choices of heat bath temperatures. 

Since also in the 2d case the temperature gradient vanishes in the thermodynamic limit, one is 
allowed to conjecture that linear response theory should reproduce the behavior of sufficiently 
large systems. 
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According to the general discussion carried on in Sec. 5.2, 
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dr, 



(183) 







where V — Ra^N^ is the system volume and J*^^^ is the x-componcnt of the total heat flux 
vector (182). Notice that simulations have to be performed for a sufficiently large size A^^, since 
the thermodynamic limit has to be taken before the infinite time limit in the above formula. 

Numerical simulations at constant energy, with R — 1/2 and periodic boundary conditions in 
both directions [69] convincingly suggest a logarithmic divergence (in time) of the correlation 
integral appearing in (183). Once more, this scaling behavior is consistent with the outcome of 
direct non-equilibrium simulations. Indeed, by assuming that Eq. (183) reproduces the correct 
size dependence if the integral is cut-off at a time t — aN^/vg [72], the temporal logarithmic 
divergence translates into an analogous divergence with the system size (see also Section 6.2). 

Let us finally mention that recent results by Shimada et al. [121] confirm the overall scenario 
and, furthermore, provide the direct confirmation that k is finite in 3d for this class of models. 



10 Conclusions 

While this review, hopefully, provides a rather complete account of the existing dynamical 
approches to heat conduction in low dimensional lattices, it certainly does not solve all open 
questions. Some of the most interesting issues requiring further investigations are briefly sum- 
marized in this concluding chapter. 

The first problem concerns heat transport at low temperatures. It is well known that many of the 
Hamiltonian models used for describing anharmonic crystals with nearest-neighbor interactions 
exhibit very slow relaxation to equilibrium below a typical energy density Cc (or temperature) 
that depends on the model and on the space dimension. For instance, the Id FPU /?-model 
shows a crossover between fast and slow relaxation at a value of the energy density e ^ 1, with 
all the parameters of the model set to unity [102]. The same holds for the Id Lennard-Jones 
potential at a close value of the energy density for e and a (see Eq. (3) ) chosen in such a 
way that the coefficients of the second and fourth order terms of the Taylor series expansion 
around the equilibrium position coincide with those of the FPU /?-model. When passing to 2d, 
for both models, the value of Cc decreases: for instance, in the Lennard-Jones 6/12 potential 
Cc ~ 0.3 [122]. This peculiar behavior in the low-temperature regime can be attributed to 
long living meta-stable states that slow-down dramatically the relaxation process . A similar 

It would be interesting to investigate the possibility of experimental tests of such a phenomenon 
in real solids, where transient effects can be usually resolved by fast spectroscopic techniques. Some 
authors have also suggested strong analogies with glassy dynamics [123], a subject that has recently be- 
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scenario can be observed in the Id rotor model described in Section 6.4. When approaching 
the two intcgrablc Umits of this model (the harmonic and the "free rotors" for small and large 
temperatures respectively) again slow relaxation mechanisms set in. An even more interesting 
situation concerns the 2d version of the rotator model, akin to the Xy-model (see e.g. [124] 
and references therein). In fact, it is characterized by the presence of the so called Kosterlitz- 
Thouless phase transition at finite temperature between a disordered high temperature phase 
and a low temperature one, where vortices condensate. From a dynamical point of view, there 
are analogies with the above mentioned examples, namely also this low temperature phase 
exhibits slow relaxation dynamics to equilibrium [119]. The possible relation with topological 
changes of the phase space due to the presence of a phase transition should be investigated. 

Upon all what we have discussed in this review, in particular in Section 5, one might expect that 
the slowing-down of relaxation processes at low temperatures should be even more relevant for 
transport properties. In fact, recent numerical investigations [69] have shown that in 2d FPU- 
like models the heat-flux correlation function seemingly exhibits the resurgence of a power-law 
divergence of the heat conductivity in the thermodynamic limit . Anyway, it is quite difficult 
to conclude only on the basis of numerical investigations if this has to be attributed to finite 
size and finite time effects (see the discussion reported in Ref. [69]). The problem could be 
better tackled in the Id FPU model, where finite-size effects are expected to be even more 
relevant below the crossover temperature. This question is closely related to other finite-size 
effects observed in these models [63] . 

Another issue that remains partially unexplored concerns the possible role of nonlinear excita- 
tions in transport. The idea that solitons may play a role in heat conduction dates back to Toda 
[109]. For instance, it has been invoked to explain the anomalous behavior of the FPU model as 
a consequence of ballistic transport due to solitons of the modified Kortcweg-deVries equation. 
Actually, such an equation can be obtained as a continuum limit of the FPU lattice model. 
Numerical experiments indicate that such solutions may persist as long living states of the 
FPU dynamics. On the other hand, upon what reported in the previous chapters, the leading 
contribution to the divergence of heat conductivity is given by the slow-relaxation properties of 
long wavelength modes. We cannot however exclude that also nonlinear excitations like solitons 
or kinks, according to the model at hand, may contribute to the divergence of k,. Recently, it 
has been proposed that transport properties should be affected also by the presence of periodic, 
spatially localized lattice waves denoted as breathers [108]. Anyway, the effect of any kind of 
nonlinear excitation is quite difficult to be detected. Even if we could assume that the the en- 
ergy flux is the sum of a phononic and a solitonic contribution, J = Jph + Jsoi how can we hope 
to distinguish the latter if the phononic part already yields anomalous behavior? In general, 
the chaotic features of the dynamics prevent the possibility of disentangling the contribution 
of nonlinear waves from that of extended modes. In this respect some better insight on the 
role of nonlinear excitations might be obtained from the analysis of integrable systems, as in 
the problem of ballistic spin transport (see [111] and references therein). Altogether, nonlinear 



come of primary interest for theoretical end experimental investigations in out-of-cquilibrium physics. 

It is worth recalling that for this class of 2d models, the heat conductivity shows a logarithmic 
divergence with the system size, see chapter 9. 
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Fig. 40. Finite-length conductivity in disordered FPU chains as from from Ref. [125], T+ = 1 • 10~^, 
r_ = 5 • 10~^. Diamonds refer to free boundary conditions while full dots refer to fixed b.c. In the 
inset we report the effective exponent defined in (106). 

excitations are one of the possible ingredients affecting the divergence of heat conductivity. 
Although it seems rather natural to conejcture that all models characterized by momentum 
conservation fall in the same universality class (with the only exception of bounded potentials), 
the numerical discrepancies among the various systems are at least suggestive of relatively 
strong finite-size corrections. If, on the other hand, we remind that there is no way to control 
the approximations implicitely contained in the self-consistent mode-coupling theory, we real- 
ize that even the problem of determining the asymptotic growth rate of heat conductivity in 
homogeneous systems may have not yet come to an end. 

A further remark concerns the combination of disorder and nonlinearity. This has been con- 
sidered sporadically in the past [32,113,117,118] in relation with the problem of heat transport 

(see also section 9.1). Wc should say that poor progress has been made in this direction during 
the last decades. In fact, few results are available and in general no clear conclusion can be 
drawn about the dependence of the thermal conductivity on the system size. Even the most 
recent results contained in a contribution by Hu et al. [125] are not convincing in this respect. 
In Fig. 40 we compare the finite-length heat conductivity k{N) for free and fixed boundary 
conditions in the same model and for the same parameter values considered in [125]. Upon 
these results, one can only conclude that the dynamical regime explored in that paper is practi- 
cally indistinguishable from the disordered harmonic case (see Section 4.2). Despite the already 
heavy numerical efforts needed to produce the data in Fig. 40, much longer time scales and 
system sizes have to be explored in order to fully appreciate the role of nonlinear terms. 

Wc conclude this section by addressing the reader to a final interesting open question about the 
study of heat conduction in structurally disordered lattices, as well as in models of amorphous 
materials [116] or quasicrystals [126]. When the crystal structure of a lattice is destroyed, the 
phononic contribution to anomalous heat transport is expected to play a much less relevant role. 
Nonetheless, peculiar transport properties are known to arise in real amorphous materials [127]. 
The most studied examples are real glasses, where viscosity exhibits a dramatic increase below 
a transition temperature specific of the material at hand. It is still unclear if phenomena like 
this should be ascribed to mechanisms other than phononic contributions. The effectiveness of 
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mode-coupling theory in describing thermodynamic-limit divergences in models of solids, as well 
as glassy dynamics in models of real glasses, indicate the extremely fascinating perspective of 
a possible unified theory of anomalous transport in condensed matter systems. In this respect, 
we should also remark that much remains to be done in both cases in order to clarify the 
reliability of the analytical estimates based on the mode-coupling approach. In particular, the 
many approximations adopted in the derivation of the scaling laws reported in section 5.3 are 
justified by quite rough arguments. A closer inspection of their validity by performing analytical 
as well as more accurate numerical calculations would be highly desirable. Furthermore, the 
problem of thermodynamic-limit divergences is not the only open question in this context. The 
other crucial facet of the problem is the temperature dependence of heat conductivity in such 
materials, where various dynamical solutions associated with the non-homogeneous structure 
of the systems should be considered as responsible of deviations from the expected classical 
laws. 
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A A rigorous definition of temperature 

In this appendix we discuss a rigorous dynamical definition of temperature. The starting point is 
the entropy S since in the /^-canonical ensemble, it plays the role of a generalized thermodynamic 
potential which allows determining (through the computation of suitable derivatives) any other 
thermodynamic observable. 

In particular, the temperature can be defined from the well known thermodynamic relation. 



where the subscript V indicates partial derivative at constant volume. 

Upon assuming that the phase space is equipped with a uniform undecomposable probability 
measure, S is given by the logarithm of the volume covered by all micro-states with energy 





n<E, 
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S{E, N, V) = In Q{E, N, = In j dV, 

H<E 



(A.2) 



where we have neglected an irrelevant multiplicative factor in front of f2 necessary only to make 
the argument of the logarithm dimensionless. 

From Eqs. (A.l), (A.2) one obtains 

H=E " " 



where the integral is the "area" of the constant energy hyper-surface Ti, = E. Upon now 
introducing a vector u such that V • u = 1 and using the divergence theorem, we obtain 



1 



r ij<7 f da 

_ •'n=E ||VH|| _ Jn=E \\s/n\\ . . 



It is easy to recognize that the above equation coincides with Eq. (6) which can thus be obtained 
by following a purely geometrical approach. More details about the derivation of this formula 
can be found in [20]. 



B Exact solution for the homogeneous harmonic chain 



In this appendix we closely follow the procedure adopted in Ref. [51] to solve the Fokker-Planck 
(58) equation for a homogeneous harmonic chain. Starting from the equilibrium solution (67), 
let us define 



^= 2A ^ • 

Prom Eq. (64), it follows that U, V and Z satisfy the equations. 



Z = -Z^ 
V=UG+ZR 
2S - VR - RV = u[GZ - ZG] 



(B.2) 
(B.3) 
(B.4) 
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where v = uJ^ j}? is the only, dimensionless, parameter that matters. In addition, U and V are 
required to be symmetric. From the pecuhar structure of the matrices R and S, it follows that 
the l.h.s. of Eq. (B.4) is a bordered matrix (i.e., its only nonvanishing elements are located on 
the external columns and rows). Accordingly, the r.h.s. must be bordered as well, i.e. in the 
bulk, Z commutes with G. The most general structure of a matrix commuting with G in the 
interior is the linear combination of a matrix Mf^ with equal elements along the diagonals (z+j 
constant) and a matrix with equal elements along the cross-diagonals ii — j constant). The 
antisymmetry requirement for Z (see Eq. (B.2)), implies that no contribution of the second 
type is present and, more precisely, that 

Zij = Hj - i) (B.5) 

with the further constraint 0(j) = — 0(— j). The quantities 0(j) are fixed by equating the border 
elements of the commutator [G, Z] (multiplied by u) with those of the l.h.s. of Eq. (B.4), 

i/0(j) = Sji - Vij = 5ji + VN,N-j+i (B.6) 

where (p^ = hy definition. Prom Eq. (B.3) and its transposed expression it follows that U 
satisfies a similar relation to that for Z, 

GU - UG = RZ + ZR. (B.7) 

Accordingly, also U commutes with G in the bulk. The different symmetry property of U with 
respect to Z implies, however, that U is constant along the cross-diagonals. It is easy to verify 
that a solution of Eq. (B.7) is given by 

rr _(<P{t+J-l) ift+J<N 

" I (j){2N iii+j>N ■ ^ ' 

In principle, this is not the only solution of Eq. (B.7), as one can add any symmetric matrix 
commuting with G; however, one can check a posteriori that the addition of any such matrix 
would eventually violate the symmetry properties of V. 

As a result of Eq. (B.8), also the matrix X can be expressed in terms of the auxiliary variables 
0(j). By replacing the Z and X solutions in the r.h.s. of Eq. (B.3), we both obtain an equation 
for the vector 0(j), 

N-l 

^ Ki^<f>{j) = 5u, (B.9) 

where K. — G + ul, and the following expression for V, 

V = S - i/U. (B.IO) 
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The problem of finding a solution for the heat transport in a homogeneous chain is accordingly 
reduced to solving Eq. (B.9) that can be written as the recursive relation 

00- + 1) = + 2)0(i) - - 1) (B.ll) 



which has to be complemented by suitable initial and final conditions. Prom the above equation, 
it follows that is the linear combination of two exponentials exp(±Q;_7) with 




e-^l+ Ju+- (B.12) 



Upon imposing the appropriate initial conditions, we finally obtain 



which completes the solution for the stationary probability distribution. 
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